Migração por quadrados mínimos

Prévia do material em texto

UNIVERSIDADE FEDERAL DO PARÁ
INSTITUTO DE GEOCIÊNCIAS
PROGRAMA DE PÓS-GRADUAÇÃO EM GEOFÍSICA
DISSERTAÇÃO DE MESTRADO
Migração por quadrados mínimos: estudo de
estratégias de formulação e precondicionamento
JOSÉ JÉSUS DA SILVA SOBRINHO
Belém � Pará
2020
JOSÉ JÉSUS DA SILVA SOBRINHO
Migração por quadrados mínimos: estudo de
estratégias de formulação e precondicionamento
Dissertação apresentada ao Programa de Pós-Graduação
em Geofísica do Instituto de Geociências da Universi-
dade Federal do Pará para obtenção do título de Mestre
em Geofísica.
Área de concentração:
Modelagem e inversão de dados geofísicos
Linha de pesquisa:
Imageamento de dados geofísicos
Orientador: Dr. Jessé Carvalho Costa
Belém � Pará
2020
Dados Internacionais de Catalogação na Publicação (CIP) de acordo com ISBD
Sistema de Bibliotecas da Universidade Federal do Pará
Gerada automaticamente pelo módulo Ficat, mediante os dados fornecidos pelo(a) autor(a)
S586m Silva Sobrinho, José Jésus da.
 Migração por quadrados mínimos: estudo de estratégias de
formulação e precondicionamento / José Jésus da Silva Sobrinho. —
2020.
 51 f. : il. color.
 Orientador(a): Prof. Dr. Jessé Carvalho Costa
 Dissertação (Mestrado) - Universidade Federal do Pará,
Instituto de Geociências, Programa de Pós-Graduação em
Geofísica, Belém, 2020.
 1. Imageamento sísmico. 2. Migração por quadrados
mínimos. 3. Precondicionamento. I. Título.
CDD 550
Powered by TCPDF (www.tcpdf.org)
Let everything happen to you
Beauty and terror
Just keep going
No feeling is �nal.
-Rainer Maria Rilke
AGRADECIMENTOS
Gostaria de primeiramente agradecer à minha família que me forneceu todo suporte para
que eu chegasse até aqui: meu pai, que me deu todo apoio quando decidir sair para
estudar; minha mãe, pelas orações e por todo carinho e amor; minha irmã caçula, que
esquenta o meu coração e com a qual eu relembro todos os dias a felicidade absoluta que
a gente só consegue encontrar nos olhos de uma criança; e a minha irmã Apoenna, que
esteve do meu lado ao longo de todo esse caminho durante os bons e maus momentos.
Isso é por vocês.
Em segundo lugar, gostaria de agradecer o apoio e amizade dos muitos amigos que �z
e/ou fortaleci durante esses anos de dedicação ao mestrado: aos parceiros de café, Luize,
Adriany, Natiê, Isadora; aos parceiros de copo, Sidônio, Cláudia, Murilo, Caique, Iury,
Luan, Joaquim e muitos outros. Dos amigos que o mestrado me trouxe, eu gostaria de
agradecer em especial ao meu orientador professor Jessé, e ao meu parceiro de programa-
ção Bruno, pois ambos ajudaram e incitaram muitas das re�exões que vieram a resultar
neste trabalho.
Agradeço também todo apoio �nanceiro que recebi durante toda minha vida acadê-
mica desde a graduação: a CAPES, pela oportunidade do intercâmbio durante a gra-
duação com o Ciência sem Fronteiras e pela bolsa de mestrado; a SBGf, pela bolsa de
iniciação cientí�ca durante a graduação; e a PETROBRAS, pelo apoio �nanceiro ao pro-
jeto que possibilitou essa pesquisa. Ainda no âmbito acadêmico, agradeço a todo apoio
do Programa de Pós-Graduação em Geofísica da UFPa.
Por último, mas mais importante que tudo, agradeço a Deus por ter me dado força e
saúde para chegar até aqui, e por ter colocado as pessoas certas nesta minha caminhada.
RESUMO
A migração reversa no tempo por mínimos quadrados (LSRTM, na sigla em inglês) é
capaz de produzir imagens sísmicas com resolução mais alta, amplitudes relativas balan-
ceadas e reduzir os artefatos causados pela aquisição. Isto é atingido através do ajuste
dos dados de re�exão observados com um modelo linearizado para o campo de onda es-
palhado. Os dados calculados podem ser obtidos por diferentes métodos. Neste trabalho,
dois diferentes modelos para o campo de onda re�etido são investigados: a aproximação de
Born e a aproximação de Kirchho�. As diferenças e vantagens dessas duas abordagens são
estudadas através de experimentos numéricos usando dados sintéticos. Adicionalmente,
investigamos um precondicionador com o objetivo de acelerar a convergência do método
de gradientes conjugados e produzir imagens com amplitudes relativas mais balanceadas.
Palavras-chaves: Migração por quadrados mínimos. LSRTM. Aproximação de
Born. Aproximação de Kirchho�. Precondicionamento.
ABSTRACT
Least-squares reverse time migration (LSRTM) produces seismic images with higher
resolution, improved relative amplitude balance, and reduced artifacts produced by acqui-
sition footprints. This is achieved by �tting the observed re�ection data with a linearized
model for the scattered wave�eld. The linearized model data might be obtained th-
rough di�erent methods. In this work, two di�erent models for the re�ected �eld are
investigated: the Born approximation and the Kirchho� approximation. The di�erences
between the approaches and the advantages of each model are studied through numerical
experiments using synthetic data sets. Additionally, we investigate a preconditioner to
accelerate the convergence of the conjugated gradient method and to produce images with
more balanced relative amplitudes.
Keywords: Least squares migration. LSRTM. Born approximation. Kirchho� ap-
proximation. Preconditioning.
LISTA DE FIGURAS
2.1 (a) Ponto de re�exão em um re�etor plano. (b) Ponto de re�exão em um
re�etor com mergulho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Ilustração esquemática do RTM em duas dimensões espaciais. Começando
com o máximo tempo registrado tmax, o campo de onda é injetado a partir
da superfície de observação em direção a subsuperfície, então é retropropa-
gado no tempo até t0. A equação da onda propaga energia descendente e
ascendente dentro do domínio computacional. A cada intervalo temporal
um corte vertical ao longo do cubo mostra um instantâneo do campo de
onda. Este campo de onda é então cross-correlacionado com o campo da
fonte, que é propagado forward no tempo, segundo a condição de imagem
utilizada. Adaptado de Etgen et al. (2009). . . . . . . . . . . . . . . . . . . 5
2.3 Diagrama esquemático da geometria para os gradientes do tempo de trân-
sito da fonte da função de Green, ∇τs, e o gradiente do receptor da função
de Green, ∇τr. Onde θ é o ângulo de abertura entre o raio incidente e o
raio re�etor. Adaptado de Yao et al. (2018). . . . . . . . . . . . . . . . . . 9
3.1 (a) Modelo de re�etividade verdadeiro. (b) Imagem do LSRTM com apro-
ximação Born depois de 20 iterações. (c) Imagem do LSRTM com aproxi-
mação Kirchho� depois de 20 iterações. (d) Primeira derivada na direção
z da imagem do LSRTM com aproximação de Born. . . . . . . . . . . . . . 19
3.2 4 primeiros shotgathers do dado observado para o modelo Marmousi2. . . . 21
3.3 (a) Wavelet broadband utilizada nas modelagens com banda de frequência
em Hz: [2.0, 10.0, 22.0, 30.0]. (b) Espectro de amplitudes da wavelet. . . . 21
3.4 (a) Modelo de velocidade verdadeiro do Marmousi2. (b) Modelo de veloci-
dade suavizado do Marmousi2. (c) Modelo de re�etividade do Marmousi2. 22
4.1 Iterações do LSRTM utilizando a aproximação de Born sem precondicio-
nador: (a) iteração 1, (b) iteração 25. . . . . . . . . . . . . . . . . . . . . . 24
4.2 Iterações do LSRTM utilizando a aproximação de Kirchho� sem precondi-
cionador: (a) iteração 1, (b) iteração 25. . . . . . . . . . . . . . . . . . . . 25
4.3 Iterações do LSRTM utilizando a aproximação de Born com precondicio-
nador Laplaciano: (a) iteração 1, (b) iteração 25. . . . . . . . . . . . . . . 26
4.4 Iterações do LSRTM utilizando a aproximação de Kirchho� com precondi-
cionador Laplaciano: (a) iteração 1, (b) iteração 25. . . . . . . . . . . . . 27
4.5 Iterações do LSRTM utilizando a aproximação de Born com precondicio-
nador iluminação: (a) iteração 1, (b) iteração 25. . . . . . . . . . . . . . . 28
4.6 Iterações do LSRTM utilizando a aproximação de Kirchho� com precondi-
cionador iluminação: (a) iteração 1, (b) iteração 25. . . . . . . . . . . . . .29
4.7 Curvas de resíduo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.8 Resultados �nais do LSRTM depois de 25 iterações: (a) aproximação de
Born e precondicionador Laplaciano, (b) aproximação de Born e precon-
dicionador iluminação, (c) aproximação de Kirchho� e precondicionador
Laplaciano, (d) aproximação de Kirchho� e precondicionador iluminação,
(e) re�etividade verdadeira do modelo Marmousi2. . . . . . . . . . . . . . . 31
4.9 Espectros de amplitude do: (a) LSRTM com aproximação de Born e pre-
condicionador Laplaciano (b) LSRTM com aproximação de Born e pre-
condicionador iluminação, (c) LSRTM com aproximação de Kirchho� e
precondicionador Laplaciano, (d) LSRTM com aproximação de Kirchho�
e precondicionador iluminação, (e) modelo de re�etividade verdadeiro do
modelo Marmousi2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
SUMÁRIO
1 INTRODUÇÃO 1
2 FUNDAMENTAÇÃO TEÓRICA 3
2.1 MIGRAÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1 Migração Reversa no Tempo . . . . . . . . . . . . . . . . . . 5
2.1.2 Migração Como um Problema Inverso . . . . . . . . . . . . 6
2.2 OPERADORES DO PROBLEMA DIRETO . . . . . . . . . . . . . . . . 7
2.2.1 Aproximação de Born . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.2 Aproximação de Kirchho� . . . . . . . . . . . . . . . . . . . . 8
2.3 MÉTODO DOS ESTADOS ADJUNTOS . . . . . . . . . . . . . . . . . . 11
2.4 PRECONDICIONANDO O SISTEMA LINEAR . . . . . . . . . . . . . . 16
3 METODOLOGIA 18
3.1 IMPLEMENTAÇÃO DO LSRTM . . . . . . . . . . . . . . . . . . . . . . 18
3.2 APROXIMAÇÃO BORN X APROXIMAÇÃO KIRCHHOFF: COMPA-
RAÇÃO EM UM MODELO DE CAMADAS HORIZONTAIS . . . . . . . 19
3.3 MODELO MARMOUSI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4 RESULTADOS 23
4.1 SEM PRECONDICIONADOR . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 PRECONDICIONADOR LAPLACIANO . . . . . . . . . . . . . . . . . . 23
4.3 PRECONDICIONADOR ILUMINAÇÃO . . . . . . . . . . . . . . . . . . 24
4.4 ESPECTROS DE AMPLITUDE E COMPARAÇÃO DOS RESULTADOS 26
5 CONCLUSÃO 33
REFERÊNCIAS 35
APÊNDICES 38
A� GRADIENTE CONJUGADO PRECONDICIONADO E REGULARI-
ZADO 39
1 INTRODUÇÃO
A migração por quadrados mínimos é reconhecida como um importante método para a
estimativa da re�etividade da terra, capaz de produzir imagens com uma qualidade maior
do que aquelas produzidas por métodos convencionais de migração (Zhang et al., 2014). A
teoria da migração por quadrados mínimos foi primeiramente derivada por Nemeth et al.
(1999) para aplicações na migração Kirchho�; posteriormente, com o avanço do poder
computacional, aplicações da migração por quadrados mínimos foram realizadas com o
método de migração reversa no tempo. Atualmente, a migração reversa no tempo é o
estado-da-arte dos métodos de imageamento, principalmente para áreas com estruturas
complexas na subsuperfície, e quando a mesma é empregada como um problema inverso
de quadrados mínimos é capaz de reduzir ruídos devido a limitações dos algoritmos de
migração em aquisições sísmicas com parametrizações insatisfatórias, atenuar artefatos
de baixo número de onda e mitigar os efeitos de iluminação em subsuperfície, resultando
em uma imagem com maior resolução (Dai and Schuster, 2013).
Na migração por quadrados mínimos a relação linearizada entre o dado calculado e
o parâmetro invertido pode ser deduzida por duas formas: a primeira, mais comumente
encontrada na literatura, utiliza a aproximação de Born - foi a abordagem utilizada por
Nemeth et al. (1999) quando primeiro formulou a migração por quadrados mínimos - e
tem como parâmetro físico invertido pertubações do modelo de impedância; a segunda
alternativa utiliza a aproximação de Kirchho� para a linearização do problema, resultando
em um parâmetro de inversão relacionado a re�etividade (Yao et al., 2018).
Além disso, na migração por quadrados mínimos as atualizações do modelo são cal-
culadas a partir da minimização da função objetivo. Portanto, um passo intermediário
mandatório na migração por quadrados mínimos é a avaliação do gradiente da função
objetivo. Entretanto, o gradiente construído a partir da operação de correlação cruzada
sofre com efeitos do espalhamento geométrico, o que resulta em amplitudes baixas para
re�etores profundos e uma taxa de convergência mais lenta (Huang et al., 2016). Xu and
Sacchi (2017) tentaram mitigar esse problema aplicando um precondicionador à função
objetivo que compensasse pela energia do campo de onda da fonte e obtiveram uma acele-
ração na convergência e uma melhor resolução espacial. Já Huang et al. (2016) utilizaram
uma aproximação da matriz Hessiana como um precondicionador para compensar pelo
efeito do espalhamento geométrico e melhorar a performance da inversão por quadrados
mínimos.
Portanto, motivado por tais considerações, este trabalho tem como objetivos princi-
pais: 1) fórmular a LSRTM com as duas aproximações lineares - Born e Kirchho� - para
assim analisar suas diferenças e possíveis vantagens do uso de uma ou outra; e 2) propor
um precondicionador que, idealmente, acelere a taxa de convergência do método de gradi-
1
2
ente conjugado utilizado e, ao mesmo tempo, entregue imagens com amplitudes relativas
dos re�etores mais próximas possíveis das amplitudes relativas do modelo de re�etividade
verdadeiro; uma vez que a obtenção de imagens com amplitudes relativas preservadas é
essencial na análise de litologias (Wang, 2017). Para tal �m, nós analisamos o efeito de
dois precondicionadores diferentes na migração por quadrados mínimos; sendo o primeiro
um operador Laplaciano convencional, uma vez que é bem conhecida sua propriedade de
atenuar ruídos de baixo número de onda; e o segundo um precondicionador para compen-
sar pela iluminação do campo de onda da fonte, sendo tal operador inspirado no inverso
assintótico da modelagem Born derivado por Op't Root et al. (2012).
Esta dissertação está estruturada da seguinte maneira: no capítulo 2 apresentamos
o referencial teórico onde analisamos as motivações que levaram ao desenvolvimento do
método RTM, a formulação da migração como um problema inverso, onde deduzimos as
duas formulações estudadas para o operador direto, demonstramos a dedução do cálculo
da função objetivo pelo método dos estados adjuntos, e �nalizamos esse capítulo apre-
sentando os precondicionadores analisados neste trabalho. No capítulo 3 tratamos da
implementação do LSRTM e analisamos as diferenças na utilização da aproximação de
Born ou Kirchho� nos resultados do LSRTM em um modelo simples de três camadas ho-
rizontais. No capítulo 4 analisamos os resultados das duas formulações do LSRTM com a
aplicação dos precondicionadores propostos em experimentos realizados no modelo Mar-
mousi2. Em seguida o capítulo 5 apresenta a discussão de nossos resultados e conclusões.
Por �m, o apêndice A� contém a dedução do método de gradiente conjugado regularizado
utilizado para a solução do problema linear.
2 FUNDAMENTAÇÃO TEÓRICA
2.1 MIGRAÇÃO
Na prospecção de hidrocarbonetos o principal método geofísico utilizado é a sísmica
de re�exão, onde o campo de onda gerado por uma fonte na superfície é re�etido em
descontinuidades geológicas na subsuperfície e então registrado em receptores localizados
na superfície. Esses dados são então processados para a obtenção de uma seção sísmica
que corresponda idealmente a uma seção com as principais interfaces do local estudado,
levando-se em consideração alguns pressupostos. No empilhamento, que é uma das etapas
do processamento, pode-se pressupor que as camadas em subsuperfície são horizontais.
Logo, para um par fonte-receptor com afastamento não-nulo, o ponto re�etor se encon-
traria diretamente abaixo do ponto médio entre a fonte e o receptor (Figura 2.1a). Então
pode-se organizar todo o conjunto de dados em famílias de ponto médio comum para a
realização da etapa de correção NMO e em seguida o empilhamento, obtendo-se,assim, a
seção empilhada que é uma aproximação de uma seção de afastamento nulo (seção zero-
o�set). Tal seção já pode ser considerada uma primeira aproximação de uma imagem da
subsuperfície.
No entanto, para o caso de camadas inclinadas o ponto re�etor não estará localizado
diretamente abaixo do ponto médio entre a fonte e o receptor (Figura 2.1b), logo a seção
empilhada resultará em uma imagem na qual os re�etores inclinados serão focalizados em
posições diferentes da posição verdadeira. Tal deslocamento na real posição de um ponto
de um re�etor inclinado é corrigido por meio de uma migração.
A migração, com a posse de um modelo de velocidade su�cientemente acurado, obje-
tiva mover re�etores inclinados para suas verdadeiras posições em subsuperfície ao mesmo
tempo que colapsa difrações, portanto aumentando a resolução espacial e fornecendo uma
imagem sísmica da subsuperfície (Yilmaz, 2001). Em outras palavras, o objetivo da mi-
gração é fazer com que a seção sísmica �nal corresponda a uma imagem com as principais
interfaces do local estudado. Nesta etapa os pontos re�etores aparentes são reposicionados
para os pontos re�etores verdadeiros em profundidade ou tempo. A denominação migra-
ção é devido a esse reposicionamento das coordenadas dos pontos re�etores aparentes para
os pontos re�etores verdadeiros.
Na realidade, a migração é mais antiga que a maioria dos métodos de processamento
sísmico. Durante a década de 1920 surgiram as primeiras implementações da migração
como um método grá�co, e antes do advento do processamento digital a migração fora
desenvolvida em diversas formulações baseadas no princípio cinemático do empilhamento
em curvas de difrações, que é o princípio por detrás da migração de Kirchho� (Gray et al.,
2001). Com o desenvolvimento do empilhamento de Ponto Médio Comum (CMP, na sigla
3
4
(a)
(b)
Figura 2.1: (a) Ponto de re�exão em um re�etor plano. (b) Ponto de re�exão em um
re�etor com mergulho.
em inglês) e a aplicação de técnicas de processamento de sinal digital em dados sísmicos
nos anos 1960, foi então possível o surgimento dos primeiros métodos de migração baseados
na equação da onda: migração como solução de diferenças �nitas de uma aproximação da
equação da onda (Claerbout and Doherty, 1972), a migração Kirchho� (Schneider, 1978)
e a migração no domínio f-k (Stolt, 1978).
Há uma variedade de estratégias de migração: pré ou pós-empilhamento; migração em
tempo ou em profundidade; migração 2-D ou 3-D. Em suma, o espectro de estratégias de
migração estende-se desde migração 2-D pós-empilhamento em tempo até migração 3-D
pré-empilhamento em profundidade. Dependendo da geologia da subsuperfície do local
estudado, qualquer estratégia entre estes extremos pode ser selecionada (Yilmaz, 2001).
A estratégia utilizada neste trabalho é a Migração Reversa no Tempo (RTM, na sigla em
inglês) que é um tipo de migração pré-empilhamento em profundidade.
5
2.1.1 Migração Reversa no Tempo
Na Migração Reversa no Tempo (Baysal et al., 1983; McMechan, 1983; Whitmore,
1983) o campo de onda observado é propagado a partir das fronteiras do domínio (a
superfície de observação) em direção ao interior da Terra utilizando a equação completa da
onda com o tempo regredindo, o que é possível devido à simetria temporal da equação da
onda. Para o caso da migração pré-empilhamento, o campo da fonte é também calculado
utilizando-se a equação da onda, mas com o tempo na direção normal (forward time),
em seguida, o campos da fonte e dos receptores são utilizados no cálculo da condição de
imagem. A Figura 2.2 ilustra com mais detalhes todo o processo.
É importante notar que o método RTM é baseado no teorema da reciprocidade do tipo
correlação (Nemeth et al., 1999; Rosa, 2010). No entanto, por motivos práticos não é pos-
sível adquirir observações do campo de onda em uma superfície fechada como idealizado
no teorema da reciprocidade. Ainda assim, mesmo com observações limitadas e discre-
tas na superfície o método RTM produz resultados satisfatórios, conseguindo focalizar
pontos re�etores ou difratores com relativa precisão dependendo do modelo de velocidade
utilizado.
Figura 2.2: Ilustração esquemática do RTM em duas dimensões espaciais. Começando
com o máximo tempo registrado tmax, o campo de onda é injetado a partir da superfície de
observação em direção a subsuperfície, então é retropropagado no tempo até t0. A equação
da onda propaga energia descendente e ascendente dentro do domínio computacional. A
cada intervalo temporal um corte vertical ao longo do cubo mostra um instantâneo do
campo de onda. Este campo de onda é então cross-correlacionado com o campo da fonte,
que é propagado forward no tempo, segundo a condição de imagem utilizada. Adaptado
de Etgen et al. (2009).
6
A versatilidade da RTM vem do fato desta ser um método de migração baseado na
solução da equação da onda completa, devido a isto ela é capaz de lidar com variações
laterais de velocidade e produzir muitos fenômenos ondulatórios (re�exões, refrações, di-
frações, múltiplas e ondas evanescentes). Além disso, na RTM a propagação da onda se dá
em todas as direções possíveis, e sua habilidade de imagear mergulhos com forte inclinação
lhe dá uma vantagem sobre outros métodos (Bednar, 2005; Etgen et al., 2009). Atual-
mente, para muitos projetos de imageamento cujo propósito seja delinear trapas sutis de
hidrocarbonetos sob um meio sobrejacente complexo a RTM é considerada indispensável.
O uso da solução da equação da onda, em vez de uma aproximação, é capaz de dar ao
RTM uma melhoria na formação da imagem quando comparado com outros métodos de
migração baseados na teoria do raio ou nas aproximações parabólicas da equação da onda.
Além disso, uma modi�cação simples na condição de imagem pode possibilitar ao RTM
gerar imagens compensadas pelo espalhamento geométrico, resultando em amplitudes re-
lativas mais �éis para os re�etores mais profundos (Xu et al., 2011; Hou and Symes, 2015;
Dafni and Symes, 2018). Discutimos este tópico na seção 4.3 pois um dos objetivos deste
trabalho é analisar possíveis precondicionadores para o LSRTM que forneçam imagens
mais próximas possíveis das verdadeiras amplitudes relativas.
A utilização da solução da onda completa na formulação da RTM ocasiona um pro-
blema particular que não é observado em outros métodos de migração, como os baseados
na teoria do raio, que é o surgimento de artefatos de baixo número de onda acima de
regiões com alto contraste de impedância. Esses artefatos surgem devido ao fato dos
campos de onda da fonte e dos receptores terem o mesmo caminho de propagação ao
longo do qual a energia correlaciona; então, a correlação cruzada desses dois campos não
produz amplitude apenas nos re�etores, mas também em todos pontos não-re�etores ao
longo de todo o percurso entre fonte e receptor, resultando em ruídos de alta amplitude
que contaminam a imagem próximo de interfaces com forte contraste de impedância (Liu
et al., 2011). Na seção 2.4 discutimos alguns precondicionadores estudados neste trabalho
que podem ser utilizados para mitigar os efeitos dos artefatos de baixo número de onda.
2.1.2 Migração Como um Problema Inverso
O objetivo �nal da Geofísica é fazer inferência sobre as propriedades físicas da sub-
superfície a partir dos dados observados (Zhdanov, 2015). Tal processo é habitualmente
denominado de inversão e é expresso matematicamente pela seguinte equação:
d = L(m), (2.1)
onde d é o dado observado, m é o modelo desconhecido e L é o operador que descreve
a relação entre o dado e o modelo da subsuperfície. No caso da inversão sísmica, há
diferentes atributos sísmicos que podem ser invertidos: tempo de trânsito, informação da
7
fase e informação da forma de onda (Schuster, 2017). Para efeitos deste trabalho, nós
estamos interessados na inversão da forma de onda dos eventos de re�exão.
A equação 2.1 trata de problemas não-lineares, que são os mais comuns em geofísica
e os mais difíceisde resolver. Devido a isto procuram-se formas mais simples de se resol-
verem problemas não-lineares através de aproximações destes por problemas lineares. No
caso da migração, Tarantola (1984) foi um dos primeiros a formulá-la como um problema
inverso linear, neste caso a equação 2.1 pode ser reescrita como:
d = Lm, (2.2)
onde d é o vetor N × 1 dos dados observados, m é o vetor M × 1 dos parâmetros do
modelo, e L é o operador linear que mapeia do espaço dos modelos para o espaço dos
dados e tem dimensão de uma matriz N × M . A equação 2.2 representa o problema
direto, isto é, ela permite o cálculo dos dados a partir dos parâmetros m do modelo;
logo L é também conhecido como operador de modelagem do problema direto. No caso
do LSRTM, o operador de modelagem nos fornece os dados de re�exão primárias. Na
próxima seção abordaremos duas maneiras de se obter tal operador: a primeira a partir
da aproximação de Born (Tarantola, 1984; Dai et al., 2012; Hou and Symes, 2016), e uma
segunda que utiliza a aproximação de Kirchho� (Yao et al., 2018; Trad, 2020).
A migração é o convencionalmente de�nida como o adjunto do operador de modelagem,
pois ela mapeia do espaço dos dados de re�exões para o espaço do modelo de pertubação
(Claerbout, 1992; Chen and Sacchi, 2017; Wu and Bai, 2018). Então, isso é representado
da seguinte forma:
m = LTd, (2.3)
onde LT representa o operador de migração que faz a projeção dos dados no espaço dos
modelos. Na seção 2.3 deduzimos o operador LT a partir do método dos estados adjun-
tos. Com os operadores de modelagem e adjunto em mãos temos todas as ferramentas
necessárias para a realização do LSRTM.
2.2 OPERADORES DO PROBLEMA DIRETO
2.2.1 Aproximação de Born
A migração por quadrados mínimos pressupõe que se possa estimar a re�etividade
da subsuperfície a partir dos registros do campo de onda re�etido considerando apenas
espalhamentos de primeira ordem em um meio de fundo suave. Para a equação da onda
acústica, temos:
1
c(x)2
∂2p (x, t;xs)
∂t2
−∇2p (x, t;xs) = s (t;xs) , (2.4)
A aproximação de Born assume que o modelo de velocidade, c(x), pode ser separado
8
em uma componente do meio de fundo suave, c0(x), e uma componente de pertubação,
δc(x), i.e.,
c(x) = c0(x) + δc(x) , (2.5)
onde |δc(x)| << c0(x). Desta maneira o problema direto é linearizado, logo eventos de
re�exões múltiplas não são considerados (Tarantola, 1984).
De forma parecida, o campo de pressão, p (x, t;xs), pode ser decomposto como
p (x, t;xs) = p0 (x, t;xs) + δp (x, t;xs) , (2.6)
no qual p0 (x, t;xs) é o campo de onda que se propaga no meio de fundo com velocidade
suave que, portanto, obedece
1
c0(x)2
∂2p0 (x, t;xs)
∂t2
−∇2p0 (x, t;xs) = s (t;xs) , (2.7)
para um pulso fonte s (t;xs) injetado na posição xs.
Substituindo as expressões 2.5 e 2.6 na equação 2.4, e mantendo somente os termos
de primeira ordem, nós obtemos a equação da onda linearizada
1
co(x)2
∂2δp (x, t;xs)
∂t2
−∇2δp (x, t;xs) =
m(x)
c2o(x)
∂2po (x, t;xs)
∂t2
, (2.8)
onde m(x) = 2δc(x)/co(x) representa a re�etividade.
Formalmente, a solução de diferenças �nitas da equação (2.8) pode ser representada
em uma notação matricial da seguinte forma
Lm = d , (2.9)
onde L representa o operador de modelagem Born, m é o modelo de re�etividade e d
representa o campo de onda espalhado de primeira ordem δp(x, t;xs), ou dado Born,
utilizado para ajustar o campo de onda re�etido observado.
2.2.2 Aproximação de Kirchho�
Além da aproximação de Born utilizada para avaliar o dado calculado, podemos utilizar
a integral de Kirchho� para calcularmos os dados de re�exão (problema direto). Nesta
seção seguimos a abordagem de Yao et al. (2018) para deduzir o operador direto utilizando
a aproximação de Kirchho�.
O campo de onda incidente em um ponto de uma interface re�etora pode ser expresso
como:
pi(x, ω) = Gs(x|xs, ω)s(ω), (2.10)
9
onde x representa a coordenada do ponto re�etor, Gs(x|xs, ω) é a função de Green da
posição da fonte para a posição x, e s(ω) é o pulso fonte.
O fenômeno de re�exão tem como característica mudar a direção do campo de onda
incidente e escalar a amplitude do mesmo pela re�etividade, rk, mantendo assim a forma
de onda do campo incidente inalterada. Logo, o campo de onda re�etido pode ser expresso
como:
pr(x, ω) = rk(x, α)Gs(x|xs, ω)s(ω), (2.11)
onde rk(x, α) é a re�etividade com dependência angular, e α é igual ao ângulo de re�exão
que é igual à metade de θ (Figura 2.3).
Uma vez que o campo de onda re�etido é conhecido na interface re�etora, o campo de
onda no receptor pode ser computado a partir da integral de Kirchho�
ps(xr, ω) =
∫
S
(
Gr(x|xr, ω)
∂pr(x, ω)
∂n
− pr(x, ω)
∂Gr(x|xr, ω)
∂n
)
dS, (2.12)
onde xr é a posição do re�etor, Gr(x|xr, ω) é a função do Green da posição do receptor
para x, S é a interface re�etora, e ∂pr
∂n
representa a derivada direcional de pr na direção
normal:
∂pr
∂n
≡ ∇pr · n (2.13)
A derivada direcional do campo re�etido, pr(x, ω), ao longo do vetor normal à interface
Figura 2.3: Diagrama esquemático da geometria para os gradientes do tempo de trânsito
da fonte da função de Green, ∇τs, e o gradiente do receptor da função de Green, ∇τr.
Onde θ é o ângulo de abertura entre o raio incidente e o raio re�etor. Adaptado de Yao
et al. (2018).
10
re�etora, n, pode ser escrita como (a partir da Eq. (2.11)):
∂pr(x, ω)
∂n
= −rk(x, α)
∂Gs(x|xs, ω)
∂n
s(ω), (2.14)
a equação (2.14) é conhecida como aproximação de Kirchho�.
O campo de onda re�etido pr(x, ω) e o campo de onda incidente se propagam em
direções opostas, logo a presença do sinal negativo. No contexto da teoria do raio, quando
o comprimento de onda típico do sinal é muito menor que o comprimento de onda típico
das pertubações no modelo de velocidade, a função de Green pode ser aproximada como
(Cerveny, 2005):
G(x, ω) ≈ A(x)e−iωτ(x|xs), (2.15)
onde A representa a amplitude e τ o tempo de trânsito. Logo, sua derivada direcional é
dada por:
∂G
∂n
≈ −iωn · ∇τG, (2.16)
onde ∇τ é o gradiente do tempo de trânsito.
Substituindo-se as equações (2.14) e (2.16) na equação (2.12):
ps(xr, ω) = iω
∫
S
n. (∇τs +∇τr)Gr(x|xr, ω)rk(x, α)Gs(x|xs, ω)s(ω)dS, (2.17)
A equação (2.17) descreve a contribuição de um único re�etor para o dado de re�exão,
ps(xr, ω). O dado de re�exão para todo o modelo pode ser obtido integrando sob a direção
normal n:
ps(xr, ω) = iω
∫ ∫
S,n
n · (∇τs +∇τr)Gr(x|xr, ω)rk(x, α)Gs(x|xs, ω)s(ω)dSdn, (2.18)
onde n é a variável de integração ao longo da direção normal a interface.
A partir da Figura 2.3, temos que para o meio isotrópico:
∇τs +∇τr =
2 cosαn
v(x)
, (2.19)
onde α é o ângulo de re�exão e v(x) é a velocidade de propagação na posição x.
Substituindo-se a equação (2.19) na equação (2.18)
ps(xr, ω) = iω
∫ ∫
S,n
2 cosα
v(x)
Gr(x|xr, ω)rk(x, α)Gs(x|xs, ω)s(ω)dSdn, (2.20)
Transformando a integral da equação (2.20) sob n e S para x:
ps(xr, ω) = iω
∫
2 cosα
v(x)
Gr(x|xr, ω)rk(x, α)Gs(x|xs, ω)s(ω)dx, (2.21)
11
Se a imagem da re�etividade empilhada (ângulo independente) for para ser recuperada,
então a equação (2.21) pode ser simpli�cada como:
ps(xr, ω) = iω
∫
Gr(x|xr, ω)Ir(x)Gs(x|xs, ω)s(ω)dx, (2.22)
onde Ir(x) é a imagem da re�etividade empilhada. De acordo com o teorema da recipro-
cidade, a equação (2.22) pode ser reescrita como
ps(xr, ω) = iω
∫
Gr(xr|x, ω)Ir(x)Gs(x|xs, ω)s(ω)dx, (2.23)
A equação (2.23) pode ser entendida como solução de um sistema de equações no
domínio do tempo:
1
v02(x)
∂2p0(x, t)
∂t2
−∇2p0(x, t) =
∂s(t)
∂t
, (2.24)
1
v02(x)
∂2ps(x, t)
∂t2
−∇2ps(x, t) = Ir(x)p0(x, t). (2.25)
O sistema de Equações (2.24) e (2.25) são muito parecidos com o sistema deduzido
na seção anterior para o cálculo do dado Born, a diferença está na injeção da fonte
do campo de onda no modelo de referência e no escalonamento da injeção da fonte do
campo re�etido. Devido a isto, é possível implementar as duas abordagens com o mesmo
modelador, fazendo poucas modi�cações de uma abordagem para a outra.
2.3 MÉTODO DOS ESTADOSADJUNTOS
O LSRTM consiste na solução do problem de quadrados mínimos linearizado
argmin
m
∥∥Lm− dobs∥∥2
2
(2.26)
onde dobs é o dado sísmico de re�exão observado. Neste trabalho esse problema inverso é
resolvido através de iterações de gradiente conjugado (GC). Portanto, para implementar
essa solução nós precisamos do gradiente da função objetivo (2.26), que corresponde a um
operador migração aplicado à diferença entre os dados observados e calculados (Tarantola,
1984; Claerbout, 1992). Nesta seção nós deduzimos uma expressão para esse gradiente
por meio do método dos estados adjuntos (Plessix, 2006).
Em muitos problemas de inversão em Geofísica estamos interessados em minimizar o
erro entre os dados observados e o dado calculado a partir de um modelo do problema es-
tudado. No caso da migração por quadrados mínimos estamos interessados em minimizar
o resíduo entre os dados observados e modelados para um conjunto de fontes e receptores,
12
então a função objetivo que tentamos otimizar tem a forma:
J (ps;m) =
1
2
∑
s
∑
t
∑
r
[
pobs (t,xr;xs)− ps (t,xr;xs|m(x))
]2
, (2.27)
onde pobs é o dado de re�exão observado e ps o dado calculado (Born) que obedece a
equação 2.8. O objetivo é minimizar J vinculado à equação de evolução do campo Born
2.8:
1
co(x)2
∂2ps (x, t;xs)
∂t2
−∇2ps (x, t;xs) =
m(x)
c2o(x)
∂2po (x, t;xs)
∂t2
. (2.28)
Determinar um extremo de J sujeito à equação de evolução acima pode ser entendido
como uma extensão do método dos multiplicadores de Lagrange (Chavent, 2010) para
determinar o extremo de uma função sujeita à restrições. Em cada instante de tempo
e em cada posição do domínio no qual ocorre a propagação as equações de evolução
devem ser satisfeitas, os multiplicadores de Lagrange devem ser especi�cados para cada
instante e coordenada (t,x); ou seja, os multiplicadores de Lagrange passam a ser campos
adjuntos de�nidos em todo o domínio de propagação durante todo o intervalo de tempo
de modelagem, portanto é um campo associado (adjunto) à equação de evolução 2.28.
Analogamente à minimização de uma função, especi�camos um funcional estendido,
a Lagrangeana, formado pela adição de J ao vínculo de�nido pela equação de evolução
ponderada pelo campo adjunto, especi�camente:
L (ps,Λ|m) =
1
2
∑
s
∑
t
∑
r
[
pobs (t,xr;xs)− ps (t,xr;xs|m)
]2
+
∑
s
∫
Ω(x)
dΩ
∫ T
0
dtΛ (t,x;xs)
[
1
co(x)2
∂2ps (x, t;xs)
∂t2
−∇2ps (x, t;xs)−
m(x)
c2o(x)
∂2po (x, t;xs)
∂t2
]
,
(2.29)
na qual Ω representa o domínio de propagação, o intervalo de registro dos campos é de 0
a T ; e Λ é o estado adjunto ao campo de pressão.
Para determinar o gradiente de J em relação a perturbações na propriedade física,
neste caso a re�etividade m(x), avaliamos variação de primeira ordem do funcional da
equação (2.29) em relação a uma con�guração estacionária do campo de onda. Conside-
rando uma perturbação na re�etividade, avaliamos o efeito de primeira ordem no funcional
em torno de uma con�guração de extrema dos campos de onda, ou seja, con�gurações que
satisfazem as equações de campo (2.28). Para avaliar a variação do funcional consideramos
13
as con�gurações de campo e perturbações do parâmetro físico:
p′s(t,x) = ps(t,x) + δps(t,x),
m′(x) = m(x) + δm(x).
(2.30)
Adicionalmente pressupomos que estas variações sejam nulas na fronteira do domínio e
que os campos de onda satisfação as mesmas soluções iniciais. Substituindo as expressões
acima na Lagrageana (2.29) e retendo apenas os termos de primeira ordem obtem-se:
δL =
∑
s
∫
Ω
dΩ(x)
{
−
∑
r
δ (x− xr)
[∫ T
0
dt
[
pobs(t,x)− ps(t,x)
]
δps(t,x)
}
+
∑
s
∫
Ω(x)
dΩ
∫ T
0
dtΛ (t,x;xs)
[
1
co(x)2
∂2δps (x, t;xs)
∂t2
−∇2δps (x, t;xs)−
δm(x)
c2o(x)
∂2po (x, t;xs)
∂t2
]
.
(2.31)
Considerando-se as parcelas de δL ponderadas por Λ(t,x),
∑
s
∫
Ω(x)
dΩ(x)
∫ T
0
dtΛ(t,x;xs)
[
1
c20(x)
∂2δps
∂t2
−∇2δps −
δm(x)
c2o(x)
∂2po
∂t2
]
. (2.32)
Integrando por partes a parcela que envolve a segunda derivada em relação ao tempo
das variações no campo de onda ps:∫ T
0
dtΛ(t,x)
∂2δps
∂t2
= Λ(t,x)
∂δps
∂t
∣∣∣∣T
0
−
∫ T
0
dt
∂Λ
∂t
∂δps
∂t
. (2.33)
Considerando a condição �nal Λ(T,x) = 0, temos:∫ T
0
dtΛ(t,x)
∂2δps
∂t2
= −
∫ T
0
dt
∂Λ
∂t
∂δps
∂t
. (2.34)
Integrando por partes novamente o lado direito da equação acima:
−
∫ T
0
dt
∂Λ
∂t
δps
∂t
=
∂Λ
∂t
δps
∣∣∣∣T
0
+
∫ T
0
dt
∂2Λ(t,x)
∂2t
δps. (2.35)
E tendo como condições �nais ∂Λ(T,x)
∂t
= 0 e Λ(T,x) = 0, temos:∫ T
0
dtΛ(t,x)
∂2δps
∂t2
=
∫ T
0
dtδps
∂2Λ(t,x)
∂2t
. (2.36)
Analisando agora a parte que envolve o laplaciano das variações do campo de onda
ps na equação (2.32), vemos que podemos utilizar a regra da cadeia para reescrever esse
14
termo como: ∫
Ω(x)
dΩΛ(t,x)∇2δps =
∫
Ω(x)
dΩ
[
∇ · (Λ∇δps)−∇Λ · ∇ps
]
. (2.37)
Aplicando-se o Teorema de Gauss no lado direito da equação acima, obtemos:∫
Ω(x)
dΩΛ(t,x)∇2δps =
∫
∂Ω
Λ∇δps · n̂ dS −
∫
Ω
dΩ∇Λ .∇δps. (2.38)
Considerando-se condições de fronteira nulas para Λ(t,x), temos:∫
Ω(x)
dΩΛ(t,x)∇2δps = −
∫
Ω(x)
dΩ∇Λ · ∇δps. (2.39)
Utilizando novamente a regra da cadeia podemos reescrever o lado direito da equação
acima como:
−
∫
Ω(x)
dΩ∇Λ · ∇δps = −
∫
Ω(x)
dΩ
[
∇ · (∇Λ δps)−∇2Λ(t,x)δps
]
. (2.40)
Pelo Teorema de Gauss, temos:
−
∫
Ω(x)
dΩ∇Λ · ∇δps = −
∫
∂Ω
dS δps ∇Λ · n̂ +
∫
Ω(x)
dΩδps∇2Λ(t,x). (2.41)
Considerando-se as condições de fronteira ∇Λ.n̂ = 0, temos:
−
∫
Ω(x)
dΩ∇Λ · ∇δps =
∫
Ω(x)
dΩδps∇2Λ(t,x). (2.42)
Logo, para as condições de fronteira Λ(t,x) = 0 e ∇Λ · n̂ = 0, temos que:∫
Ω(x)
dΩΛ(t,x)∇2δps =
∫
Ω(x)
dΩδps∇2Λ(t,x). (2.43)
Substituindo os resultados de (2.36) e (2.43) na equação de variação de primeira ordem
da Lagrangeana (2.31), podemos obter a expressão para a variação da função objetivo em
decorrência de variações no modelo m em qualquer região do domínio Ω:
δL(ps,Λ|m) =
∑
s
∫
Ω(x)
dΩ
∫ T
0
dtδps
[
1
c20(x)
∂2Λ(t,x)
∂t2
−∇2Λ(t,x)
−
∑
r
(pobs(t,x)− ps(t,x))δ(x− xr)
]
−
∑
s
∫
Ω
dΩ
∫ T
0
dtδm
Λ(t,x)
c20(x)
∂2p0(t,x)
∂t2
.
(2.44)
15
Como o objetivo é determinar as variações do funcional devidas exclusivamente às
variações da propriedade física, neste caso pertubações na impedância ou re�etividade
representadas por m, a expressão acima não pode depender de δps. Consequentemente
devemos exigir que os coe�cientes de δps se anulem, o que resulta na seguinte equação
para o campo adjunto:
1
c20(x)
∂2Λ(t,x)
∂t2
−∇2Λ(t,x) =
∑
r
(
pobs(t,x)− ps(t,x)
)
δ(x− xr). (2.45)
Sujeito às condições �nais:
Λ(T,x;xs)
∣∣
Ω
= 0 . (2.46)
∂Λ(T,x;xs)
∂t
∣∣
Ω
= 0 . (2.47)
e a condições de fronteira homogêneas em Λ e ∇Λ.n, de�nidas em consequência das
condições de fronteira do problema de modelagem e recobrindo toda a fronteira do domínio
de propagação, ∂Ω. Portanto, o estado adjunto é solução da equação acústica escalar em
que os termos fontes correspondem aos resíduos entre o campo observado e o campo
modelado.
Determinando o estado adjunto Λ, a variação da Lagrangeana pode ser imediatamente
obtida,
δL(ps,Λ|m) = −
∑
s
∫
Ω
dΩ
∫ T
0
dtδm
Λ(t,x)
c20(x)
∂2p0(t,x)
∂t2
. (2.48)
Devido a forma como foi de�nida a Lagrangeana (2.29), a variação de L devido a
pertubações no modelo m equivale a variação da função objetivo J :
δJ = δL = −
∑
s
∫
Ω
dΩ
∫ T
0
dtδm
Λ(t,x)
c20(x)
∂2p0(t,x)
∂t2
. (2.49)
Do cálculo, temos que a variação em J devido a pertubações e m(x) é:
δJ = ∇mJ · δm(x) . (2.50)
Comparando a equação (2.49) com a (2.50) e reconhecendo a integral no domínio Ω
como um produto escalar, temos que:
∇mJ = −
∑
s
∫ T
0
dtΛ(t,x)
1
c20(x)
∂2p0(t,x)
∂t2
. (2.51)
onde ∇mJ é o gradiente da função objetivo na direção da pertubação do modelo m(x).
Analisando-se a equação (2.51) percebe-se que à parte de uma multiplicação pela ve-
16
locidade do meio de fundo e de uma derivada parcial temporal de segunda ordem, o
gradiente da função objetivo é similar a uma imagem migrada e sua fórmula é cinemati-
camente similar ao princípio do imageamento (Claerbout, 1985).
2.4 PRECONDICIONANDOO SISTEMA LINEAR
A convergência do método dos gradiente conjugados pode ser comprometida por um
gradiente desbalanceado, marcas de aquisição e artefatos de baixo número de onda pre-
sentes no gradiente da função objetivo. Além disso, como apontado por Claerbout (2014),
métodos iterativos, como o das direções conjugadas, podem ser muitas vezes acelerados
por uma mudança de variáveis, onde a mudança de variável mais simples é chamada de
solução teste. Então, com objetivo de reduzir os artefatos provenientes de marcas da
aquisição e de acelerar a convergência do método do gradiente conjugado, nós precondici-
onamos o problema linear com uma mudança de variáveis. Formalmente, nós escrevemos
a solução como:
m = Pa, (2.52)
onde m é a solução que procuramos, as colunas da matriz P são formatos que preferimos
e os coe�cientes em a são coe�cientes desconhecidos para selecionar as quantidades dos
formatos prediletos. As variáveis a são comumente chamadas de variáveis precondiciona-
das.
Idealmente, o operador P deve ser projetado para �ltrar artefatos de baixo número
de onda, acelerar a taxa de convergência das iterações do GC e melhorar a qualidade do
modelo invertido. Para �ltrar artefatos de baixo número de onda o operador P não pode
ser positivo de�nido e para garantir a convergência do GC nós precisamos adicionar um
termo de regularização para a, resultando no seguinte sistema linear regularizado:(
LP
µI
)
a =
(
dobs
0
)
, (2.53)
onde µ é o parâmetro de regularização, tal parâmetro é utilizado para tornar a solução do
sistema linear mais estável; e I é a matriz identidade. O apêndice A� apresenta a dedução
do algoritmo utilizado na solução desse problema.
Nós investigamos o efeito de dois tipos de precondicionadores para o LSRTM. O pri-
meiro precondicionador é um operador Laplaciano 2D, o qual possui a propriedade de
�ltrar as componentes de baixo número de onda (Guitton et al., 2006). Neste caso, a
equação (2.52) assume a forma:
m =
[
∇2
]
a . (2.54)
Nós propomos a segundo precondicionador que foi inspirado pelo inverso assintótico
do operador da modelagem Born derivado por Op't Root et al. (2012). Este operador
17
além de aplicar um operador Laplaciano também compensa pela iluminação da fonte.
Neste caso P assume a forma:
m =
[
c2(x)∇2
I(x)
]
a , (2.55)
onde I é a iluminação do campo de onda da fonte (Kaelin and Guitton, 2006). Este
precondicionador quando aplicado ao gradiente da função objetivo (equação (2.51)) difere
da condição de imagem deduzida por Op't Root et al. (2012) por esta última ter como
fonte do problema adjunto a integral do campo observado nos receptores, mas ambas
fórmulas multiplicam pelo quadrado da velocidade e dividem pela iluminação do campo
da fonte.
3 METODOLOGIA
3.1 IMPLEMENTAÇÃO DO LSRTM
Para a realização do LSRTM é necessária a implementação de duas operações: migra-
ção (cálculo do gradiente da função objetivo) e modelagem direta (cálculo do dado). Em
ambas operações o modelo de velocidade do meio de fundo (�gura 3.4b) é utilizado, e tal
modelo é uma das entradas do LSRTM e equivale ao termo c0(x) da equação 2.5. Esse
modelo se mantém constante ao longo das iterações do LSRTM e a solução do sistema
linear é o modelo de pertubação δc(x).
Para a realização da modelagem direta e adjunta a equação da onda acústica (equação
2.4) foi discretizada no tempo com diferenças �nitas com um precisão de segunda ordem,
e o laplaciano foi avaliado utilizando também o método de diferenças �nitas com uma
precisão de sexta ordem. A condição de borda utilizada para a atenuação do campo nas
extremidades do domínio foi a C-PML (Pasalic and McGarry, 2010), que devido a sua
e�ciência permite a utilização de uma borda com poucos pontos, facilitando, desta forma,
uma economia de memória.
Como demonstrado no capítulo anterior, a diferença de implementação da modelagem
direta utilizando a aproximação de Born e Kirchho� está apenas na fonte. Para facili-
tar o desenvolvimento e manutenção do código nós o construímos com uma abordagem
inspirada no paradigma de programação orientada a objeto. Os serviços necessários -
como injeção da fonte, avaliação da condição de imagem, avanço no tempo, algoritmo de
otimização - foram programados da forma mais independente possível, o que possibilitou
a reutilização de muitas subrotinas e uma melhor organização.
As modelagens direta e adjunta apresentam uma independência com relação aos ti-
ros, isto é, os dados gerados devido a uma fonte na posição xs = x0 não têm nenhuma
relação com os dados com uma segunda fonte na posição xs = x0 + dx. Isto nos permitiu
utilizar as subrotinas de MPI (Message Pass Interface) para realizar a modelagem dos
tiros em paralelo, diminuindo, desta forma, o tempo total de modelagem de acordo com
a quantidade de núcleos disponíveis.
Com todos os serviços acimas mencionados implementados de forma independente, o
último passo para a realização do LSRTM foi a escolha do método de otimização para
a resolução do sistema linear. Para este trabalho nós escolhemos o gradiente conjugado
precondicionado e regularizado, cuja dedução encontra-se no Apêndice A�.
Todos os códigos deste trabalho foram implementados em FORTRAN 2003 e paraleli-
zados utilizando MPI 3.0, com exceção dos plots dos resultados que foram feitos utilizando
a biblioteca Matplotlib do Python (Hunter, 2007).
18
19
3.2 APROXIMAÇÃO BORN X APROXIMAÇÃO KIRCHHOFF: COMPARAÇÃO EM
UM MODELO DE CAMADAS HORIZONTAIS
Como vimos nas seções anteriores, duas aproximações podem ser utilizadas para a
modelagem direta do dado para o LSRTM: aproximação de Born ou aproximação de
Kirchho�. Elas possuem signi�cados físicos diferentes. A aproximação de Born é baseada
em pertubações na impedância, enquanto que a aproximação de Kirchho� é deduzida
usando a re�etividade (Yao et al., 2018).
0 200 400 600 800 1000
Distancia (m)
0
200
400
600
800
1000
Pr
of
un
di
da
de
 (m
)
(a)
0 200 400 600 800 1000
Distancia (m)
0
200
400
600
800
1000
Pr
of
un
di
da
de
 (m
)
 
 
(b)
0 200 400 600 800 1000
Distancia (m)
0
200
400
600
800
1000
Pr
of
un
di
da
de
 (m
)
 
 
(c)
0 200 400 600 800 1000
Distancia (m)
0
200
400
600
800
1000
Pr
of
un
di
da
de
 (m
)
 
 
(d)
Figura 3.1: (a) Modelo de re�etividade verdadeiro. (b) Imagem do LSRTM com aproxi-
mação Born depois de 20 iterações. (c) Imagem do LSRTM com aproximação Kirchho�
depois de 20 iterações. (d) Primeira derivada na direção z da imagem do LSRTM com
aproximação de Born.
20
Para compararmos o efeito da in�uência da utilização da aproximação de Kirchho�
ou de Born no resultado �nal do LSRTM, realizamos um experimento em um modelo
simples de três camadas, onde a segunda camada possui uma velocidade v = 2200m/s, e
a primeira e terceira camadas possuem uma velocidade v = 2000m/s. O modelo possui um
grid de 200x200 com dz = dx = 5m, e 9 tiros foram distribuídos igualmente na superfície
a partir de x = 50m até x = 950m. O dado foi gerado utilizando-se uma wavelet Ricker
de 30 Hz. A velocidade de migração é uma versão suavizada do modelo de velocidade
verdadeiro.
As Figuras 3.1b e 3.1c mostram o resultados do LSRTM depois de 20 iterações
utilizando-se a aproximação Born e a aproximação Kirchho�, respectivamente. Neste ex-
perimento utilizamos a iluminação da fonte para precondicionar o problema linear (Kaelin
and Guitton, 2006). Comparando-se os resultados obtidos, percebe-se que a aproxima-
ção de Kirchho� produz uma imagem simétrica (fase zero) do re�etor, enquanto que a
aproximação de Born produz uma imagem anti-simétrica do re�etor. Isso implica que a
posição do re�etor é indicado pelo pico na imagem da aproximação de Kirchho� (seta
na Figura3.1c), enquanto que, utilizando a aproximação de Born, a posição do re�etor é
indicada pela zona de transição do zero na imagem (seta na Figura 3.1b).
Se a primeira derivada ao longo da direção z é aplicada na imagem do LSRTM com
aproximação Born (Figura 3.1b),o resultado (Figura 3.1d) é muito similar a imagem
do LSRTM com aproximação de Kirchho�. Yao et al. (2018) provou que isso é devido
ao fato da imagem da aproximação de Kirchho� ser a primeira derivada da imagem da
aproximação Born sob a pressuposição de independência angular. Entretanto, essa simples
conversão, que implica na derivada ao longo da profundidade, só é válida para um modelo
de camadas horizontais.
3.3 MODELO MARMOUSI
O dado de entrada é de uma porção do modelo Marmousi2 (Martin et al., 2002)
e foi avaliado com um modelador diferente do implementado nas modelagens diretas e
adjuntas do LSRTM. Nós procuramos evitar, desta forma, o que alguns autores, como
Schuster (2017), de�nem como o crime da inversão: a utilização do mesmo modelador
para produzir o dado observado e o dado calculado. O período de gravação do dado foi
de 3.5 s com uma amostragem de 0.004 s. A �gura 3.2 mostra as quatro primeiras seções
de tiro-comum do dado sintético, com a aplicação de um silenciamento para a retirada
da onda direta, então é esperado que alguns eventos, como múltiplas, gerem artefatos que
não serão explicados pelo dado calculado, uma vez que o operador direto é um operador
linear que prediz apenas re�exões primárias.
O modelo de velocidade original foi modi�cado para uma tamanho de 9 km x 3.5 km
com um intervalo de grid de 10m (�gura 3.4a). Nos experimentos 240 tiros foram igual-
21
Figura 3.2: 4 primeiros shotgathers do dado observado para o modelo Marmousi2.
(a) (b)
Figura 3.3: (a) Wavelet broadband utilizada nas modelagens com banda de frequência em
Hz: [2.0, 10.0, 22.0, 30.0]. (b) Espectro de amplitudes da wavelet.
mente distribuídos na superfície, a posição da primeira e última fonte são, respectivamente,
3000m e 8975m, com uma distância entre cada tiro de 25m. Para cada família de tiro
comum, o o�set mínimo é 200m e o máximo 2575m, com 96 receptores igualmente dis-
postos com um espaçamento de 25m entre si. A �gura 3.3 mostra a wavelet utilizada
na modelagem e o seu respectivo espectro. A injeção da fonte pontual e os registros
dos receptores pontuais foram implementados usando uma função de banda limitada que
fornece medidas acuradas em posições arbitrárias do grid (Hicks, 2002).
22
(a)
(b)
(c)
Figura 3.4: (a) Modelo de velocidade verdadeiro do Marmousi2. (b) Modelo de velocidade
suavizado do Marmousi2. (c) Modelo de re�etividade do Marmousi2.
4 RESULTADOS
4.1 SEM PRECONDICIONADOR
As Figuras 4.1 e 4.2 mostram o resultado do LSRTM usando, respectivamente, a apro-
ximação de Born e a aproximação de Kirchho�, ambas sem precondicionador. É possível
notar a presença de artefatos de baixo número de onda, e também que a aproximação de
Kirchho� se mostra mais sensível a esse tipo de ruído. É possível também perceber que
a presença saturada de artefato de baixo número de onda prejudica o melhoramento da
imagem pelo LSRTM, pois a maior parte da melhora da imagem é na redução do mesmo
como pode ser observado nas regiões em destaque nas �guras 4.1b e 4.2b.
Como visto no capítulo anterior, quando o modelo de camadas horizontais foi ana-
lisado, é possível notar uma diferença de fase nas imagens dependendo da aproximação
utilizada para o operador direto. Tal diferença é mais perceptível comparando o resultado
da primeira iteração (Figuras 4.1a e 4.2a), onde é possível veri�car que a aproximação de
Kirchho� é mais suscetível a artefatos de baixo número de onda quando não utiliza-se um
precondicionador.
4.2 PRECONDICIONADOR LAPLACIANO
As �guras 4.3 e 4.4 mostram os resultados do LSRTM usando, respectivamente, a
aproximação de Born e a aproximação de Kirchho� com o precondicionador Laplaciano.
Comparando esses resultados com os da seção anterior é possível notar que o precondicio-
nador Laplaciano é capaz de diminuir com sucesso os artefatos de baixo número de onda.
Além disso, observando o fundo do mar nas �guras 4.3b e 4.4b, em torno do crossline
= 6000m, é possível perceber a diferença de fase devido a utilização das duas diferentes
aproximações lineares para o operador direto, para o caso da aproximação de Kirchho�
(�gura 4.4b) o re�etor se localiza sob o valor máximo de re�etividade, enquanto que para
a aproximação de Born (�gura 4.3b) o re�etor se localiza na zona de transição do valor
da re�etividade.
Apesar da melhora considerável devido apenas à utilização do precondicionador Lapla-
ciano, como pode ser observado nas regiões em destaque nas �guras 4.3b e 4.4b, quando
comparamos os resultados das �guras 4.3 e 4.4 com a re�etividade verdadeira do modelo
Marmousi2 (�gura 3.4c) percebemos que as amplitudes relativas desses resultados não
coincidem com do modelo verdadeiro.
23
24
(a)
(b)
Figura 4.1: Iterações do LSRTM utilizando a aproximação de Born sem precondicionador:
(a) iteração 1, (b) iteração 25.
4.3 PRECONDICIONADOR ILUMINAÇÃO
As �guras 4.5 e 4.6 mostram os resultados do LSRTM usando, respectivamente, a
aproximação de Born e a aproximação de Kirchho� com o precondicionador iluminação.
Como vimos na Equação 2.55 proposta para esse precondicionador, há a aplicação de
um operador Laplaciano que é capaz de reduzir consideravelmente os artefatos de baixo
número de onda resultantes da correlação cruzada utilizada da condição de imagem. En-
quanto que a correção pela iluminação da fonte possibilita a recuperação de imagens com
amplitudes relativas mais �éis ao modelo de re�etividade verdadeiro (�gura 3.4c), isto é
melhor observado nas regiões mais rasas do modelo em torno do crossline = 3000m.
25
(a)
(b)
Figura 4.2: Iterações do LSRTM utilizando a aproximação de Kirchho� sem precondici-
onador: (a) iteração 1, (b) iteração 25.
Quando analisamos o efeito da aproximação linear utilizada no operador direto perce-
bemos que devido a diferença de fase entre as duas aproximações os resultados apresentam
um pequeno deslocamento vertical na localização dos re�etores, sendo que o operador que
utiliza a aproximação de Kirchho� é capaz de entregar resultados com resoluções ligei-
ramente melhores, devido ao fato de que com essa aproximação a localização do re�etor
se dá no máximo valor da re�etividade (como descrito na seção 3.2). Além disso, com-
parando a primeira iteração (�guras 4.5a e 4.6a) com a última iteração (�guras 4.5b e
4.6b) podemos observar o aumento da resolução das imagens com os re�etores mais bem
de�nidos, como pode ser observado nas regiões em destaque nessas �guras.
26
(a)
(b)
Figura 4.3: Iterações do LSRTM utilizando a aproximação de Born com precondicionador
Laplaciano: (a) iteração 1, (b) iteração 25.
4.4 ESPECTROS DE AMPLITUDE E COMPARAÇÃO DOS RESULTADOS
Nesta seção analisaremos como a convergência do método de gradientes conjugados
é in�uenciada pelos precondicionadores e pela aproximação linear utilizada (Born e Kir-
chho�), compararemos os resultados �nais do LSRTM para todos os precondicionadores
utilizados e analisaremos os espectros de amplitudes deste resultados para identi�carmos
a banda de frequência recuperada.
A �gura 4.7 mostra o resíduo da função objetivo que é um parâmetro comumente utili-
zado para analisarmos a convergência do método de gradientes conjugados. Analisando-a
percebe-se que o LSRTM utilizando a aproximação de Kirchho� sem precondicionador
27
(a)
(b)
Figura 4.4: Iterações do LSRTM utilizando a aproximação de Kirchho� com precondici-
onador Laplaciano: (a) iteração 1, (b) iteração 25.
foi o caso que apresentou a pior convergência, uma vez que o valor do resíduo diminui
vagarosamente ao longo das iterações. Tal fato está relacionado à grande quantidade de
ruídos de baixo número de onda presentes nos resultados desse experimento (�gura 4.2)
resultantes da ausência de um precondicionador. É importante notar que há muito ruído
de baixo número de onda para o caso da aproximação de Born sem precondicionador tam-
bém, no entanto, para este caso o LSRTM consegue reduzir consideravelmente tal ruído,
ajustando melhor o dado observado - o que é evidenciadopelo resíduo na �gura 4.7.
O LSRTM com aproximação de Kirchho� e precondicionador Laplaciano foi o expe-
rimento que apresentou a maior redução do resíduo (�gura 4.7). Entretanto, quando
28
(a)
(b)
Figura 4.5: Iterações do LSRTM utilizando a aproximação de Born com precondicionador
iluminação: (a) iteração 1, (b) iteração 25.
comparamos tal experimento com o resultado �nal do LSRTM utilizando a aproximação
de Kirchho� com o precondicionador iluminação, �ca evidente que apenas a análise da
evolução do resíduo não é um parâmetro ideal para de�nirmos uma con�guração ótima
para o LSRTM. Pois este último experimento (�gura 4.8d) foi o que apresentou o melhor
resultado para as amplitudes relativas dos re�etores ao longo de todo o modelo - fato evi-
denciado quando comparamos os resultados da �gura 4.8 com a re�etividade verdadeira
(�gura 4.8e) - e uma redução dos artefatos e das marcas de aquisição que resultaram numa
imagem com maior resolução.
Na �gura 4.8 está a síntese de todo este trabalho, uma vez que a mesma mostra o re-
29
(a)
(b)
Figura 4.6: Iterações do LSRTM utilizando a aproximação de Kirchho� com precondici-
onador iluminação: (a) iteração 1, (b) iteração 25.
sultado �nal do LSRTM depois de 25 iterações para os dois precondicionadores e as duas
aproximações lineares estudadas. Nela é possível observar que ambas aproximações, Born
e Kirchho�, entregam um melhor resultado quando utilizado o precondicionador ilumina-
ção (�guras 4.8b e 4.8d), pois as amplitudes relativas, principalmente dos re�etores mais
profundos, condizem mais com as amplitudes relativas do modelo de re�etividade verda-
deiro. Quando analisamos os resultados com o precondicionador Laplaciano, no entanto,
percebemos que a ausência de um termo para compensar pela iluminação nesse precon-
dicionador faz com que os re�etores mais rasos apresentem amplitudes relativas muito
mais elevadas do que a realidade. Logo, daí percebe-se a importância da compensação
30
Figura 4.7: Curvas de resíduo.
pela iluminação da fonte para a obtenção de resultados mais próximos dos modelos reais e
que, idealmente, possibilitem a análise de amplitudes para, assim, inferirmos informações
sobre a composição das camadas.
A comparação dos espectros de amplitude dos resultados �nais do LSRTM (�gura 4.9)
com o espectro de amplitude do modelo de re�etividade verdadeira (�gura 4.9e) aponta
para a hipótese de que a maior resolução obtida com a utilização do precondicionador
iluminação é devido a sua maior capacidade de reduzir artefatos de baixo número de onda,
uma vez que para a banda de frequência de 0− 25 rad/km os espectros de amplitude do
LSRTM com esse precondicionador assemelham-se mais com o do modelo de re�etividade
verdadeira. É possível notar que o LSRTM utilizando o precondicionador Laplaciano foi o
que obteve a maior banda de frequência, o que geralmente implicaria em maior resolução;
no entanto, o maior comprometimento das baixas frequências com artefatos resultou em
uma perda de resolução.
31
(a) (b)
(c) (d)
(e)
Figura 4.8: Resultados �nais do LSRTM depois de 25 iterações: (a) aproximação de Born e
precondicionador Laplaciano, (b) aproximação de Born e precondicionador iluminação, (c)
aproximação de Kirchho� e precondicionador Laplaciano, (d) aproximação de Kirchho�
e precondicionador iluminação, (e) re�etividade verdadeira do modelo Marmousi2.
32
(a) (b)
(c) (d)
(e)
Figura 4.9: Espectros de amplitude do: (a) LSRTM com aproximação de Born e pre-
condicionador Laplaciano (b) LSRTM com aproximação de Born e precondicionador ilu-
minação, (c) LSRTM com aproximação de Kirchho� e precondicionador Laplaciano, (d)
LSRTM com aproximação de Kirchho� e precondicionador iluminação, (e) modelo de
re�etividade verdadeiro do modelo Marmousi2.
5 CONCLUSÃO
Neste trabalho analisamos o efeito do precondicionamento do LSRTM no domínio do
dado com objetivo de propormos um precondicionador que entregasse idealmente imagens
com amplitudes relativas mais próximas o possível do modelo verdadeiro. Para isso,
comparamos os resultados de dois precondicionadores, Laplaciano e Iluminação, com os
resultados do problema sem precondicionador.
O precondicionador Laplaciano mostrou-se capaz de reduzir substancialmente os arte-
fatos de baixo número de onda presentes no problema sem precondicionador. No entanto,
as amplitudes relativas, principalmente dos re�etores mais rasos, não se assemelharam
com as do modelo verdadeiro. Contudo, o precondicionador que propusemos, que leva em
consideração a iluminação da fonte, apresentou os melhores resultados com imagens com
amplitudes relativas mais próximas do modelo verdadeiro de re�etividade, e uma redução
considerável dos artefatos e marcas de aquisição, resultando em uma imagem com maior
resolução. Além disso, os espectros de amplitude dos resultados �nais do LSRTM evi-
denciam que a melhor resolução entregue pelo precondicionador iluminação vêm de sua
maior e�cácia em reduzir os artefatos de baixo número de onda.
Outro objetivo deste trabalho foi analisar o efeito de duas das aproximações lineares
que podem ser utilizadas na dedução do operador direto: aproximação de Born e Kir-
chho�. Do ponto de vista da implementação, essa comparação foi facilitada por o código
ser construído inspirado no paradigma de programação orientada a objeto, o que nos pos-
sibilitou uma melhor separação dos serviços independentes comuns a ambas deduções da
modelagem direta e dos diferentes precondicionadores estudados.
A implementação do LSRTM utilizando a aproximação de Kirchho� sem precondici-
onador mostrou-se mais suscetível a ruídos de baixo número de onda, o que resultou em
toda melhora proporcionada pela migração por quadrados mínimos ser observada apenas
na redução parcial desses artefatos, em vez da recuperação de uma imagem com amplitu-
des relativas mais balanceadas e maior resolução. No entanto, esta mesma implementação
com a aproximação de Kirchho� quando realizada com o precondicionador iluminação pro-
duziu o melhor resultado, uma vez que o termo Laplaciano presente no precondicionador
iluminação obteve êxito na retirada dos artefatos e o termo de compensão da iluminação
resultou em amplitudes relativas mais próximas as do modelo de re�etividade verdadeiro.
Além disso, é preciso enfatizar a diferença de fase nos resultados produzidos pelas duas
formulações estudadas para o operador direto. A formulação utilizando a aproximação de
Kirchho� produz imagens nas quais o valor máximo de re�etividade encontra-se na posição
do re�etor, enquanto que nas imagens produzidas pela outra formulação a localização
do re�etor encontra-se no valor nulo da amplitude (para re�etores horizontais). Esta
é a provável razão para a melhor resolução da imagem produzida pelo LSRTM com a
33
34
aproximação de Kirchho�, uma vez que re�etores muito próximos podem se sobrepor mais
facilmente quando utilizando a aproximação de Born dependendo da banda de frequência
da wavelet da fonte.
REFERÊNCIAS
Baysal, E., D. D. Koslo�, and J. W. Sherwood, 1983, Reverse time migration: Geophysics,
48, 1514�1524.
Bednar, J. B., 2005, A brief history of seismic migration: Geophysics, 70, 3MJ�20MJ.
Cerveny, V., 2005, Seismic ray theory: Cambridge university press.
Chavent, G., 2010, Nonlinear least squares for inverse problems: theoretical foundations
and step-by-step guide for applications: Springer Science & Business Media.
Chen, K., and M. D. Sacchi, 2017, Elastic least-squares reverse time migration via line-
arized elastic full-waveform inversion with pseudo-hessian preconditioningelastic lsrtm:
Geophysics, 82, S341�S358.
Claerbout, J., 2014, Geophysical image estimation by example: Lulu. com.
Claerbout, J. F., 1985, Imaging the earth's interior: Blackwell scienti�c publications
Oxford.
���, 1992, Earth soundings analysis: Processing versus inversion: Blackwell Scienti�c
Publications London.
Claerbout, J. F., and S. M. Doherty, 1972, Downward continuation ofmoveout-corrected
seismograms: Geophysics, 37, 741�768.
Dafni, R., and W. W. Symes, 2018, Accelerated acoustic least-squares migration, in SEG
Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists,
4291�4295.
Dai, W., P. Fowler, and G. T. Schuster, 2012, Multi-source least-squares reverse time
migration: Geophysical Prospecting, 60, 681�695.
Dai, W., and G. T. Schuster, 2013, Plane-wave least-squares reverse-time migration: Ge-
ophysics, 78, S165�S177.
Etgen, J., S. H. Gray, and Y. Zhang, 2009, An overview of depth imaging in exploration
geophysics: Geophysics, 74, WCA5�WCA17.
Gray, S. H., J. Etgen, J. Dellinger, and D. Whitmore, 2001, Seismic migration problems
and solutions: Geophysics, 66, 1622�1640.
Guitton, A., B. Kaelin, and B. Biondi, 2006, Least-squares attenuation of reverse-time-
migration artifacts: Geophysics, 72, S19�S23.
Hicks, G. J., 2002, Arbitrary source and receiver positioning in �nite-di�erence schemes
using kaiser windowed sinc functions: Geophysics, 67, 156�165.
Hou, J., and W. W. Symes, 2015, An approximate inverse to the extended born modeling
operatoran approximate inverse operator: Geophysics, 80, R331�R349.
���, 2016, Accelerating extended least-squares migration with weighted conjugate gra-
dient iteration: Geophysics, 81, S165�S179.
Huang, Y., R. Nammour, and W. Symes, 2016, Flexibly preconditioned extended least-
squares migration in shot-record domain: Geophysics, 81, S299�S315.
35
36
Hunter, J. D., 2007, Matplotlib: A 2d graphics environment: Computing in Science &
Engineering, 9, 90�95.
Kaelin, B., and A. Guitton, 2006, Imaging condition for reverse time migration, in SEG
Technical Program Expanded Abstracts 2006: Society of Exploration Geophysicists,
2594�2598.
Liu, F., G. Zhang, S. A. Morton, and J. P. Leveille, 2011, An e�ective imaging condition
for reverse-time migration using wave�eld decomposition: Geophysics, 76, S29�S39.
Martin, G. S., K. J. Marfurt, and S. Larsen, 2002, Marmousi-2: An updated model for the
investigation of avo in structurally complex areas, in SEG Technical Program Expanded
Abstracts 2002: Society of Exploration Geophysicists, 1979�1982.
McMechan, G. A., 1983, Migration by extrapolation of time-dependent boundary values:
Geophysical prospecting, 31, 413�420.
Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration of incomplete
re�ection data: Geophysics, 64, 208�221.
Op't Root, T. J., C. C. Stolk, and V. Maarten, 2012, Linearized inverse scattering based
on seismic reverse time migration: Journal de mathématiques pures et appliquées, 98,
211�238.
Pasalic, D., and R. McGarry, 2010, Convolutional perfectly matched layer for isotropic and
anisotropic acoustic wave equations, in SEG Technical Program Expanded Abstracts
2010: Society of Exploration Geophysicists, 2925�2929.
Plessix, R.-E., 2006, A review of the adjoint-state method for computing the gradient
of a functional with geophysical applications: Geophysical Journal International, 167,
495�503.
Rosa, A. L. R., 2010, Análise do sinal sísmico: SBGF-Sociedade Brasileira de Geofísica,
Rio de Janeiro, Brasil.
Schneider, W. A., 1978, Integral formulation for migration in two and three dimensions:
Geophysics, 43, 49�76.
Schuster, G. T., 2017, Seismic inversion: Society of Exploration Geophysicists.
Stolt, R. H., 1978, Migration by fourier transform: Geophysics, 43, 23�48.
Tarantola, A., 1984, Linearized inversion of seismic re�ection data: Geophysical prospec-
ting, 32, 998�1015.
Trad, D., 2020, Assumptions and goals for least squares migration: Geophysical Prospec-
ting, 68, 802�814.
Wang, X., 2017, Relative �delity processing of seismic data: Methods and applications:
John Wiley & Sons.
Whitmore, N. D., 1983, Iterative depth migration by backward time propagation, in SEG
Technical Program Expanded Abstracts 1983: Society of Exploration Geophysicists,
382�385.
Wu, J., and M. Bai, 2018, Incoherent dictionary learning for reducing crosstalk noise in
37
least-squares reverse time migration: Computers & Geosciences, 114, 11�21.
Xu, L., and M. D. Sacchi, 2017, Preconditioned acoustic least-squares two-way wave-
equation migration with exact adjoint operator: Geophysics, 83, S1�S13.
Xu, S., Y. Zhang, and B. Tang, 2011, 3d angle gathers from reverse time migration:
Geophysics, 76, S77�S92.
Yao, G., N. V. da Silva, and D. Wu, 2018, Forward modelling formulas for least-squares
reverse-time migration: Exploration Geophysics, 49, 506�518.
Yilmaz, Ö., 2001, Seismic data analysis: Processing, inversion, and interpretation of
seismic data: Society of exploration geophysicists.
Zhang, Y., L. Duan, and Y. Xie, 2014, A stable and practical implementation of least-
squares reverse time migration: Geophysics, 80, V23�V31.
Zhdanov, M. S., 2015, Inverse theory and applications in geophysics: Elsevier.
APÊNDICES
A� GRADIENTE CONJUGADO PRECONDICIONADO E
REGULARIZADO
Nesta seção derivaremos o esquema de gradiente conjugado precondicionado e regula-
rizado utilizado na solução da LSM. A equação 2.53 pode ser reescrita como:(
LP
µI
)T (
LP
µI
)
a =
(
LP
µI
)T (
dobs
0
)
, (A�1)
De�nindo-se L̂ = LP, temos que:(
L̂T L̂ + µ2IT I
)
a = L̂Tdobs (A�2)
O esquema de gradiente conjugado baseado na equação A�2 é apresentado no algoritmo
a seguir:
Algorithm 1 Gradiente conjugado regularizado e precondicionado
a0 = 0; r0 = dobs
g0 = L̂
T r0; p0 = g0
for n = 0, 1, 2... do
qn = L̂pn
αn =
∥∥gn∥∥2∥∥qn∥∥2+µ2∥∥pn∥∥2
an+1 = an + αnpn
rn+1 = rn − αnqn
gn+1 = L̂
T rn+1
if η‖rn+1‖2 < ‖r0‖2 then
exit
end if
βn =
∥∥gn+1∥∥2∥∥gn∥∥2
pn+1 = gn+1 + βnpn
end for
m = Pa
No algoritmo acima, gn é o n-ésimo vetor gradiente regularizado, pn é o n-ésimo vetor
de direção conjugada, rn é o n-ésimo vetor de resíduo do dado (para a iteração zero é o
dado em si), qn é o n-ésimo vetor do dado modelado, αn e βn são coe�cientes de peso,
e η é um limite de acurácia pré-determinado. Este esquema cálcula implicitamente as
operações matriz-vetor, ou seja, as matrizes nunca são formadas explicitamente; em vez
disso, elas são substituídas por subrotinas de modelagem adjunta e direta que têm como
saída apenas o vetor resultande (Nemeth et al., 1999).
39
40
Como discutido na seção 2.4, o precondicionador é introduzido no esquema de gra-
diente conjugado através de uma mudança de variáveis. Isso quer dizer que para cada
iteração há um pouco mais de trabalho: em vez dos operadores L e LT , nós aplicamos, res-
pectivamente, LP e PTLT . Nosso objetivo é que uma escolha adequada de preconcionador
diminua o número de iterações necessárias para a convergência do esquema de gradiente
conjugado (Claerbout, 2014). Felizmente, o custo computacional extra da aplicação do
precondicionador P não é muito quando comparado com o de L.