Estudo de formas implementacionais da equação da onda imagem para remigração em profundidade

Estudo de formas implementacionais da equação da onda imagem para remigração em profundidade

Fernando Perin Munerato^I; Jörg Schleicher^II; Amélia Novais^III

^IDepartamento de Engenharia de Petróleo, Faculdade de Engenharia Mecânica, UNICAMP, Rua Mendeleiev s/n, Bloco C, 3° andar, Cidade Universitária, Barão Geraldo, Caixa Postal 6122, 13083-970 Campinas, SP, Brasil. Tel. [phone]; Fax [phone] -E-mail [e-mail]

^IIDepartamento de Matemática Aplicada, Instituto de Matemática, Estatística e Computação Científica, UNICAMP, Rua Sérgio Buarque de Holanda, 651, Cidade Universitária, Barão Geraldo, Caixa Postal 6065, Sala 129, 13083-859 Campinas, SP, Brasil. Tel. [phone]; Fax: [phone] -E-mail: [e-mail]

^IIIDepartamento de Matemática Aplicada, Instituto de Matemática, Estatística e Computação Científica, UNICAMP, Rua Sérgio Buarque de Holanda, 651, Cidade Universitária, Barão Geraldo, Caixa Postal 6065, Sala 140, 13083-859 Campinas, SP, Brasil. Tel. [phone]; Fax: [phone] -E-mail: [e-mail]

RESUMO

Neste trabalho estudamos teórica e numericamente novas formas implementacionais da equação da onda imagem para remigração em profundidade. Trata-se de uma equação diferencial parcial de segunda ordem semelhante à equação da onda acústica. Determinamos, além da consistência e da estabilidade, as condições de dispersão e dissipação de seis esquemas de diferenças finitas para esta equação em diversas formas obtidas por mudanças de variáveis. Estas condições não podem ser satisfeitas simultaneamente, ie., um resultado preciso não é de fácil obtenção. Testes numéricos confirmam os resultados teóricos da estabilidade para três esquemas, mas falham nos outros. Os esquemas avançados na variável da velocidade e avançados ou atrasados em profundidade foram os que tiveram os resultados mais proveitosos, pois estes esquemas satisfazem as previsões teóricas quanto à sua estabilidade e podem ser utilizados para realizar a propagação da onda imagem, sendo o primeiro para velocidades crescentes e o segundo para velocidades decrescentes.

Palavras-chave: onda imagem, remigração, diferenças finitas, análise numérica.

ABSTRACT

ln this work, we study theoretically and numerically new implementational forms of the image-wave equation for depth remigration This is a partial second-order differential equation similar to the acoustic wave equation. We determine, additionally to consistency and stability, the conditions for dispersion and dissipation of six finite-difference schemes for this equation in different forms, obtained by a change of variables. These conditions cannot be simultaneously satisfied, i.e., a precise result is not easily obtained. Numerical tests confirm the theoretical stability results for three of the investigated schemes but fail for the others. Of the tested schemes, the ones forward in the velocity variable and forward or backward in depth were the most useful ones in the sense that they satisfy the theoretical predictions about their stability and can be utilized to realize the image-wave propagation, the former for increasing velocities and the lalter for decreasing velocities.

Keywords: image-wave, remigration, finite differences, numerical analysis.

INTRODUÇÃO

Este trabalho aborda a pergunta de como resolver a equação da onda imagem para o problema de remigração em profundidade. Este problema provém da prospecção sísmica, na qual oobjetivo é a construção de uma imagem das camadas geológicas no subsolo a partir de registros de ondas que propagaram por este meio. Esta construção de uma imagem do subsolo é conhecida na sísmica como "migração". Quando é necessária a atualização da imagem do subsolo, obtida anteriormente mediante o uso de um modelo incorreto da velocidade de propagação de ondas no subsolo, correspondente a um modelo melhor, fala-se de "remigração" (Hubral et ai, 1996)

A necessidade de melhorar a imagem obtida é freqüente na sísmica, uma vez que o modelo correto da velocidade de migração não é conhecido e tem que ser determinado durante o processo da construção da melhor imagem possível. Portanto, convém variar o modelo de velocidade e construir as correspondentes imagens. Quando é feita tal variação da velocidade, observa-se que as novas imagens dos refletores (fronteiras entre as camadas geológicas) comportam-se de forma similar à propagação de uma onda. Em meios homogêneos, esta "propagação" de imagens é descrita por uma equação diferencial parcial semelhante à equação da onda (Fomel, 1994) que, consequentemente, recebe o nome de "equação da onda imagem para remigração" (Hubral et ai, 1996)

Existem equações da onda imagem para a remigração em tempo e em profundidade. Autilidade da equação da onda imagem em tempo já foi comprovada na prática (Jaya et ai, 1999). Apesar da teoria ser desenvolvida para meios homogêneos, as equações podem ser aplicadas em meios não homogêneos. Resultados recentes para a análise de velocidade, usando a remigração em tempo, também mostraram-se bastante promissores (Fomel, 2003; Novais et ai, 2008; Schleicher et ai, 2008). Neste trabalho, tratamos da equação em profundidade.

O objetivo concreto deste trabalho é a análise de métodos numéricos para solucionar a equação da onda imagem, com a finalidade de contribuir com a análise de velocidade de migração. Como foi mostrado em Fomel (2003) e Schleicher et al (2004), o uso da solução numérica das equações da onda imagem em tempo e em profundidade possui um potencial neste sentido. Para a remigração em profundidade, este potencial édifícil de ser explorado, principalmente por causa das dificuldades numéricas, especialmente a dispersão numérica dos esquemas de diferenças finitas estáveis. Aqui continuamos os estudos desta equação que foi iniciado em Schleicher et al (2003, 2004), onde foram testados esquemas explícitos do método de diferenças finitas para diversas aproximações da derivada mista assim como algumas mudanças de variáveis para a velocidade de migração. Em todos os casos estudados anteriormente não foi encontrado um método numérico que faça a remigração em profundidade através da equação da onda imagem sem que ocorresse um forte efeito de dispersão numérica, mesmo que a estabilidade e a consistência do método estejam garantidas

O efeito de dispersão ocorre quando ondas de diferentes freqüências propagam com diferentes velocidades distorcendo, desta forma, o formato do pulso da onda. Quando as velocidades das ondas de diferentes freqüências realizadas pelo esquema numérico são diferentes da solução exata do problema, fala-se de dispersão numérica. Este efeito, intrínseco de métodos de diferenças finitas, aqui é tão forte que influência a resolução nos testes numéricos quando usamos o método de diferenças finitas para simular a propagação da onda imagem. Para cada um dos esquemas investigados em Schleicher et al (2003, 2004) observou-se que ocomprimento do pulso não se mantinha igual ao comprimento obtido pela migração direta dos dados para velocidades correspondentes. Esta desvantagem do método pode ser observada quando comparamos os resultados da técnica de remigração com imagens migradas correspondentes, obtidas por migração direta dos dados sísmicos originais.

Motivados pelos resultados obtidos em Schleicher et al (2003, 2004), estudamos outros esquemas de diferenças finitas e, principalmente, outras formas desta equação. Além disso, estendemos nossa investigação para dispersão e dissipação com o intuito de minimizar estes efeitos indesejados que aparecem na remigração.

Equação da onda imagem

A equação da onda imagem é representada pela equação diferencial

onde p (x, y, v) representa o campo de onda, x e z representam as variáveis espaciais e v a velocidade de migração, supostamente constante. Nesta equação, a variável de propagação, no entanto, não é o tempo como no caso de ondas físicas convencionais, mas a velocidade de migração, v. A condição inicial para o problema é dado por um campo de ondas migrado, po. As condições de contorno são, em conseqüência da falta de dados fora do domínio aser migrado, tomadas nulas.

Observamos que em comparação com a equação da onda, ao invés da derivada segunda com respeito à variável de propagação, a equação da onda imagem contém uma derivada mista que envolve uma coordenada espacial. Por isso, o campo de onda p descrito pela Eq. (1) possui um comportamento cinemático diferente das ondas convencionais (Fomel, 1994; Hubral et al., 1996).

Mudanças de variáveis

Assim como foram testados outros esquemas de Diferenças Finitas (DF), também estudamos a possibilidade de uma mudança de variável na equação diferencial da onda imagem, com a finalidade de melhorar as propriedades de estabilidade e dispersão. As transformações mais promissoras são α=v^q (q ∈ ) ou β = ln (v/vo), onde v_ofixo representa o valor da velocidade de migração para a migração inicial. Para esses casos, chegamos nas seguintes equações diferenciais

No segundo caso, tomamos β = ln (v/vo) de forma que β seja adimensional. Observa-se que a transformação β = In(v/vo) elimina o fator que depende da variável de propagação. Infelizmente, uma transformação correspondente não é possível para a variável z, uma vez que esta introduziria outros fatores dependentes de z para o termo p_zz.

DIFERENÇAS FINITAS

Para a implementação de DF, consideramos uma malha, cujos pontos são representados da forma

,

sendo x_o, z_oe v_ovalores iniciais, Δx, Δz e Δv incrementos em x, y e v, respectivamente, e m, n e l números naturais, Além disso, p¹_m,ndenota o campo de imagem no ponto (x_m, Z_n, v_l ) da malha.

Para as derivadas segundas espaciais, usamos aproximações de quarta ordem (Strikwerda, 1989; Thomas, 1995), i.e.,

e uma expressão correspondente para p_zz.Observa-se que a escolha de esquemas de outras ordens nas derivadas espaciais não alteram os resultados principais deste trabalho, Por outro lado, a escolha de um operador de DF para aderivada mista é muito importante. No que segue, consideramos seis aproximações diferentes que podem ser implementadas como esquemas explícitos, a primeira é avançada em vez, que já foi estudada em Schleicher et al. (2003). As outras cinco são: avançada em v e atrasada em z, avançada em v e centrada em z, atrasada em v e avançada em z, atrasada em v e z e atrasada em v e centrada em z,

Usando essas aproximações para as derivadas na equação da onda imagem, Eq, (1), encontramos os seguintes esquemas de DF

• Avançado em

v

e

z:

• Avançado em

v

e atrasado em

z:

• Avançado em

v

e centrado em

z:

• Atrasado em

v

e avançado em

z:

onde: θ = 6_vlΔz/(z_n/v).

• Atrasado em

v

e

z:

• Atrasado em

v

e centrado em

z:

Como podemos observar, estes esquemas de DF podem ser considerados implícitos na velocidade. Porém, para acelerar os testes numéricos, tratamo-los como esquemas explícitos. Isto significa que precisamos de condições de contorno adicionais para iniciar os cálculos. Todas estas condições adicionais tomamos nulas.

Diferenças finitas para mudanças de variáveis

Para os casos de mudança de variável descrito pelas Eqs (2) e (3) utilizamos os correspondentes esquemas de DF. Assim, considerando uma malha como descrito acima e α_l = α_o + l Δα e β₁ = β₀+ l Δβ temos os seguintes esquemas avançados em α e z e em β e z

Caso α =v^q:

Caso β = ln (v/v₀ )

Nestes esquemas podemos observar que, para o cálculo do nível l + 1 também iremos utilizar o nível l e o próprio nível 1+1. Concluímos então que podemos resolver as equações alteradas da mesma forma que resolvemos sem fazer tais mudanças. Os outros esquemas são obtidos correspondentemente.

ANÁLISE NUMÉRICA

Para julgar se um determinado esquema de DF pode ser usado para realizações numéricas da propagação da onda imagem, é necessário investigar a consistência e a estabilidade (Teorema de Equivalência de Lax (Strikwerda, 1989; Thomas, 1995)), assim como a sua dispersão e dissipação. Para o esquema avançado em vez, a consistência, a estabilidade, a dispersão e a dissipação foram discutidas e interpretadas em Schleicher et al. (2003, 2004). Portanto, discutimos aqui os demais esquemas, utilizando os resultados já apresentados.

Consistência

Em Schleicher et al. (2003) estudamos em detalhe a consistência do esquema de DF avançado em vez. Verificamos que o esquema é consistente com a equação diferencial da onda imagem. De forma similar, podemos verificar facilmente que o mesmo é verdade para todos os cinco outros esquemas dados pelas Eqs (5) até (9) assim como para os esquemas tendo feito as mudanças de variáveis, Eqs (10) e(11)

Estabilidade

Baseado no critério de estabilidade de von Neumann (Strikwerda, 1989; Thomas, 1995), podemos verificar a estabilidade dos esquemas estudados neste trabalho. Para isso, substituímos nas equações discretizadas, Eqs (4)-(11), a seguinte expressão

onde k_xe k_zsão as freqüências espaciais (i.e., números de ondas) em x e z,respectivamente, e ξ é chamado de fator de amplificação. Um esquema de DF é estável se o seu fator de amplificação satisfaz a condição:

Quando este critério é satisfeito sob certas condições o esquema de DF é dito condicionalmente estável. Para o caso de a condição não ser satisfeita sob hipótese alguma, o esquema de DF é instável.

De forma análoga ao que foi desenvolvido em Schleicher et al. (2003) para o esquema de DF avançado em vez, aplicamos este critério de estabilidade para todos os outros esquemas. Nota-se que nesta análise, os esquemas são tratados como explícito correspondentemente à forma que foram implementados. Uma vez que os cálculos são bastante similares, resumimos aqui somente os resultados, i.e., as condições de estabilidade para cada um dos esquemas dados pelas Eqs (4) até (11)

Avançado em vez: Para efeito de comparação reescrevemos a condição de estabilidade para o caso avançado em v e z, a saber

onde

Observamos que esta condição de estabilidade dificilmente será satisfeita para um dado migrado com conteúdo de número de onda arbitrário. Se o dado contém números de onda muito baixos na direção vertical (de modo que para estes números de onda k_zΔz 1), a restrição (14) é muito severa, porque o intervalo para valores permitidos de Λ fica muito pequeno. Porém, geralmente dados realistas têm um conteúdo de número de onda limitado, permitindo assim que a condição de estabilidade (14) seja satisfeita

Especificamente, para uma malha uniforme (Δx = Δz), e sob certas hipóteses sobre oconteúdo de números de onda vertical (k_z) e horizontal (k_x) presentes nos dados a serem migrados, a condição de estabilidade se reduz a

onde v = para k_x≈ k_ze v = para kx k_z v_min é a velocidade mínima de propagação da onda imagem e z_max é a profundidade máxima do modelo. O limite inferior foi escolhido de tal forma que a estabilidade seja garantida para todo modelo, ou seja, qualquer valor de Δv tomado neste intervalo irá satisfazer a condição de estabilidade para todo z e v,

A Eq. (16) mostra que o passo será restrito e positivo, ou seja, incremento deverá estar dentro deste intervalo de estabilidade e o esquema só poderá ser utilizado para aumentar a velocidade de migração.

Observa-se que a condição k_x k_zé realista em meios com camadas predominantemente horizontais, situação geológica comum. Por outro lado, k_x≈ k_z pode se esperar em meios com camadas de mergulhos em torno de 45°. Observa-se ainda que em meios com camadas predominantemente verticais, onde k_x

k_z a condição de estabilidade dada pela Eq. (14) dificilmente será satisfeita, uma vez que o lado esquerdo estará perto de zero, Felizmente, esta é uma situação geologicamente rara.

Avançado emv e atrasado em z: Para este esquema, Eq, (5), a condição de estabilidade é dada por

onde Λ é dado pela Eq. (15). Assim, para as mesmas hipóteses nos números de ondas e condições nos incrementos espaciais, temos que

A Eq. (18) mostra que, além do passo ser restrito, ele tem que ser negativo. Em outras palavras, a velocidade sempre terá que diminuir, isto é, este esquema só pode ser utilizado para decrementar a velocidade de migração.

Avançado emv e centrado em z: De forma análoga, chegamos na seguinte condição de estabilidade para o esquema dado pela Eq. (6)

,

ou seja, não existe Λ real que satisfaça a expressão acima, pois o termo entre parênteses é sempre positivo. Portanto, o esquema considerado é instável, ou seja, sem utilidade prática.

Atrasado emv e avançado em z: Acondição de estabilidade do esquema dado pela Eq. (7) é

,

com Λ dado pela Eq. (15) a primeira condição, portanto, é sempre satisfeita se Δv > 0. A segunda fornece para uma malha uniforme

,

onde v = para k_x≈ k_ze v = para kx k_z v_min .Poderíamos pensar que este esquema permite tanto aumentar quanto diminuir a velocidade de migração. Mas isso não éverdade, porque a condição

,

contradiz a condição de consistência, pela qual Δv deve tender a 0 independentemente de Δz.

Atrasado em ve z :A condição de estabilidade do esquema dado pela Eq. (8) é

ou

a primeira condição é sempre satisfeita se Δv < 0 e a segunda fornece para uma malha uniforme

,

onde v = para k_x≈ k_ze v = para kx k_z v_min . Como no caso anterior, esta condição contradiz a condição de consistência, Logo, este esquema não pode ser usado para aumentar a velocidade de migração.

Atrasado em v e centrado em z: Aplicando o critério de estabilidade na Eq. (9), obtemos

.

Esta expressão é verdadeira para qualquer valor de Λ. Portanto, este esquema é incondicionalmente estável, i.e., permite, teoricamente, aumentar ou diminuir a velocidade de migração usando incrementos quaisquer

Estabilidade para mudanças de variáveis

Proce dendo de forma similar ao descrito acima, analisamos as condições de estabilidade para os esquemas dados pelas Eqs (10) e(11) correspondente às equações modificadas dadas pelas Eqs (2) e(3). Em tais análises, foi observada a semelhança entre as condições de estabilidade com a mudança de variável e sem tais mudanças. Para a variável α =v^>q , substitui-se o fator 4/3 nos dois termos da Eq. (15) por 4/3q. Para a variável β=ln (v/v_o), o fator 4/3 continua omesmo, adiferença agora éque o termo VI desaparece.

Vejamos as condições finais para os dois casos.

Caso α =v^q

onde v = para k_x≈ k_ze v = para kx k_z e α_min é a velocidade mínima de propagação da onda imagem tendo feito amudança de variável α = v^q. O incremento Δα deve satisfazer a condição (19) para que o esquema dado pela Eq. (10) seja estável para todos os valores de α e z envolvidos.

Caso β=ln (v/v_o)

onde v = para k_x≈ k_ze v = para kx k_z .

Como podemos ver nas Eqs (19) e (20), as condições de estabilidade para as equações da onda imagem tendo feito as mudanças de variáveis, Eqs (10) e (11), são muito parecidas com aquela encontrada sem ter feito tais mudanças de variáveis. O mesmo acontece para os outros esquemas nas novas variáveis. Desta forma, concluímos que as mudanças de variáveis não contribuem a um melhor desempenho numérico referente à estabilidade dos esquemas.

Custo computacional

A diferença mais importante entre os esquemas da equação da onda imagem usando as diferentes variáveis está no seu custo computacional, conforme a análise descrita abaixo, quando consideramos Δx = Δz e k_x = k_z.

Sejam os incrementos na variável que representa a velocidade escolhidos tal que

i.e., os incrementos são escolhidos de acordo com as condições de estabilidade (16), (19) e (20).

Para esta análise do custo computacional, consideramos as seguintes mudanças de variáveis

Assim, da Eq. (21), concluímos que

Agora, consideramos n_α, n_γ, n_β e n_vos números de passos necessários para a propagação de v_min até v_maxem cada um dos quatro casos, respectivamente, isto é.

, ,

e ,

onde α_max , α_min, β_max, β_min, γ_max e γ_min são os respectivos valores de α β e γ correspondentes às velocidades v_maxe v_min, tendo feito as mudanças de variáveis.

Logo, podemos relacionar o número de passos para um esquema com uma mudança de variável com o correspondente número do esquema original.

, e

.

Portanto, concluímos que o número de passos n_vserá sempre menor que n_αe maior que n_β e que n_γserá menor que n_vsomente se v_max> 9v_min (que na prática é raramente necessário). Esta relação é ilustrada graficamente na Figura 1.

Em conseqüência das análises feitas e da Figura 1, para todos os casos, as mudanças de variáveis serviram somente para mostrar-nos que é possível ter uma economia no tempo computacional. Neste sentido, a maior economia seria uma implementação na variável β=ln (v/v_o) (ver Fig 1).

Dispersão e dissipação numérica

Nesta seção, estudamos a principal fonte de erros numéricos do método de DF, mesmo que o esquema seja estável. Trata-se da dispersão e dissipação numérica. Estamos interessados em saber se as mudanças de variáveis melhoram essas propriedades de dispersão e dissipação.

A seguir, discutimos em detalhe o esquema avançado em α e z, Eq. (10). A discussão para os outros esquemas (Eqs (5)(9)) e a variável β (Eq. (11)) não é explicitada, uma vez que é muito similar.

É de interesse que um pulso de onda imagem p(x,z,α)deslocado para (x₁,z₁) na velocidade α₁ tenha as mesmas propriedades que possuía na sua posição original (x_o, z_o) na velocidade α_o, ou seja, esperamos que p(x₁,z₁, α₁) = p(x_o, z_o, α_o). Assim, no domínio de Fourier, conforme detalhado em Strikwerda (1989), queremos que o esquema de DF real ize um passo de propagação por Δα satisfazendo

ondeé o campo de onda imagem no domínio de Fourier, k² = k_x² +k_z²e a é a velocidade de propagação da onda imagem. Portanto, queremos que um esquema de DF ideal da equação da onda imagem multiplique a transformada de Fourier do campo de onda acada passo por e^-ikaΔw. Mas isso quase nunca acontece, pois sempre há perda de energia na propagação da onda ou as velocidades de ondas de diferentes freqüências são distintas, quando propagadas pelo esquema de DF.

O que pretendemos então é determinar quão boa é a aproximação para o termo e^-ikaΔw que o esquema de DF realiza. Para isso, consideramos um passo no esquema de DF no domínio das freqüências espaciais k_xe k_z.De acordo com a Eq. (12), sabemos que o esquema, ao realizar o cálculo de ⁿ⁺¹a partir de ⁿ multiplica o campo acada passo pelo fator de amplificação, ξ, ou seja, em símbolos.

Comparamos a Eq. (25) com a Eq. (24) para enfatizarmos suas similaridades. Para tal, escrevemos ~ como

onde c é a velocidade de fase realizada pelo esquema. É importante observar que c é uma função de Δx e Δz bem como k_x e k_z. Observamos, pela comparação das Eqs (24) e (25), que um esquema ideal seria aquele que realizasse a propagação com |ξ|= 1 e c =a.

Portanto, se c(Δx , Δz, Δα, k_x, k_z) for igual a a para todos k_xe k_zentão as ondas propagariam com a velocidade correta, mas isto não acontece para nenhum esquema de DF para equações diferenciais hiperbólicas exceto em casos triviais. Precisamos então encontrar a função c(Δx , Δz, Δα, k_x, k_z) para que possamos escolher os incrementosΔx , Δz e Δα de tal maneira que c aproxime melhor averdadeira velocidade de propagação a da equação da onda imagem, assim reduzindo a dispersão numérica.

Observamos ainda que o fator ideal e^-ikaΔw tem amplitude 1. O fato de que |ξ| ≠ 1 implica que as amplitudes variam. Se |ξ| > 1, o esquema é instável. Se |ξ| < 1, temos perda de amplitude ao longo da propagação, ou seja, dissipação. Gostaríamos então de escolher os incrementos de tal forma que |ξ| se aproxime de 1 pela esquerda.

Caso α =v^q (avançado em α e z): correspondentemente ao que foi visto para o esquema avançado em v e z, o fator de amplificação ξ do esquema avançado em a ez é dado por

onde Λα= com Λ dado pela Eq. (15), substituindo Δv por Δα.

O argumento da função exponencial é dado pelo arco-tangente da parte imaginária dividido pela parte real, isto é.

Para encontrarmos a velocidade de fasec(Δx , Δz, Δα, k_x, k_z) , dividimos a Eq. (27) por kΔα .Para melhor entender esta função, desenvolvemos o resultado em séries de Taylor de primeira ordem em Δα e de segunda ordem nos incrementos espaciais Δx eΔz. Assim, encontramos

onde k² = k_x²+ k_z². Observa-se que, nesta ordem de aproximação, c não depende de Δx. Isto se deve ao fato de que o esquema escolhido é de quarta ordem em x. Observa-se também que o efeito da dispersão é maior para os números de ondas mais altos, ou seja, valores elevados de k, e para maiores profundidades z.

Nesta aproximação para c, retemos três termos. O primeiro não depende de Δz e Δα. Ele representa a velocidade de fase realizada pela própria equação da onda imagem. Observa-se que este termo depende de k_xe k_z, isto é, a própria equação diferencial causa dispersão na propagação da onda imagem. Este efeito corresponde ao conhecido estiramento do pulso na migração (Tygel et al., 1994). Apesar de indesejado, ele não pode ser reduzido pela escolha dos valores de Δz e Δα.

Uma vez que o segundo e o terceiro termo da Eq. (28) têm sinais diferentes, é possível reduzir o seu efeito mediante cancelamento das suas contribuições. Para isso, devemos escolher Δα e Δz de tal forma que

que depende das freqüências espaciais k_xek_z, da profundidade e do valor atual da velocidade. Mesmo sob certas hipóteses sobre o conteúdo de números de onda nesses dados, diminuir este efeito não é trivial e não será possível para todas as profundidades e velocidades. A mesma conclusão pode ser tirada a partir da correspondente análise para os outros esquemas e a variável β.

Para entender a dissipação do esquema, escrevemos o módulo do fator de amplificação como um decaimento exponencial, i.e.,

Para encontrar uma expressão para o coeficiente de amortecimento b do esquema, desenvolvemos o logaritmo da amplitude do fator de amplificação ξ do esquema avançado em α e z em série de Taylor

.

Desta forma, determinamos a seguinte expressão para o coeficiente de amortecimento b:

Na Eq. (30), observa-se que todos os termos retidos dependem dos incrementos, ou seja, a equação da onda imagem da Eq. (2) não possui dissipação intrínseca. Em outras palavras, uma solução perfeita deveria preservar automaticamente as amplitudes. Além disso, observamos que b não depende somente do número de ondas k_xe k_z, mas também de z_n e α_l, i.e., a perda de amplitude é diferente em profundidades diferentes e velocidades de migração diferentes

Da mesma forma, a dissipação numérica pode ser reduzida se os incrementos podem ser escolhidos tal que o primeiro e o segundo termo da série de Taylor da Eq. (30) se cancelem. Isto conduz a

De fato, sabemos da análise da estabilidade que a condição dada pela Eq. (31) deve ser corrigida por um fator de 3/4 por causa dos termos mais elevados das séries de Taylor. Isso é porque a condição de estabilidade de von Neumann significa nada mais do que a condição de que o esquema não pode induzir crescimento exponencial. Isto equivale a dizer que o fator de amortecimento b, definido na Eq. (30), não deve tornar-se negativo. Isso é garantido se as condições de estabilidade do esquema são satisfeitas.

Análises correspondentes foram feitas para todos os outros esquemas estáveis, com e sem mudanças de variáveis. Tais análises tiveram resultados similares aos descritos acima para o esquema avançado em α e z. Em outras palavras, as mudanças de variáveis não conduzem a melhores propriedades dos esquemas neste aspecto em comparação ao esquema avançado em v e z estudado em Schleicher et al. (2004) A melhor escolha para Δα para reduzir a dispersão e a dissipação são diferentes. Isto significa que ambos efeitos numéricos não podem ser reduzidos ao mesmo tempo por uma escolha inteligente de Δα e Δz.

Além disso, a dispersão e a dissipação dependem da profundidade e da velocidade de migração. Portanto, é impossível a escolha de Δα e Δz de tal maneira que a dispersão e a dissipação sejam reduzidas para todo passo na velocidade e toda profundidade. A única maneira de se conseguir uma melhoria total de ambos os efeitos numéricos é escolher valores de Δα e Δz tão pequenos quanto os recursos computacionais permitirem. Um experimento numérico que ilustra esta propriedade, que é equivalente a demonstrar a convergência numérica do esquema, foi realizado para o esquema avançado em v e z em Schleicher et al. (2004) Este procedimento, apesar de melhorar a precisão da aproximação (diminuindo o efeito de dispersão), aumenta drasticamente o tempo computacional, pois a malha precisa ser muito mais refinada para conseguir uma boa redução da dispersão e dissipação.

RESULTADOS NUMÉRICOS

Refletor plano

Para tornar os estudos numéricos mais simples, utilizamos um modelo simples de um refletor plano horizontal em 550 m de profundidade com dados registrados na região entre -2000 m e 2000 m na direção x e entre 10 m e 1600 m na direção z, comΔx = Δz = 10 m. Um exemplo para a propagação da imagem de um refletor com o esquema avançado em v e zcomΔv > 0 já foi apresentado em Schleicher et al. (2004) assim como uma análise dos efeitos da propagação quando refinamos o incremento espacial Δz e outros tipos de refletores (inclinados e em meios inomogêneos). Aqui discutimos os resultados numéricos dos outros esquemas investigados com e sem mudanças de variáveis.

Esquema avançado em v e atrasado em z:

Este esquema é estável somente com um passo de velocidade Δv negativo. Para testá-lo, migramos os dados provenientes do experimento numérico com uma velocidade incorreta de 4000 m/s. A partir desta, remigramos a imagem do refletor plano até este atingir sua profundidade correta de 550 m com velocidade 3000 m/s (que é de conhecimento prévio já que o modelo é sintético). A Figura 2 (esquerda) ilustra a remigração, mostrando as imagens inicial, intermediária (v = 3500 m/s) e final do refletor. Escolhemos um passo na velocidade que satisfaça a condição de estabilidade dada pela Eq. (18), Δv = -1m/s Testes com passos maiores, que não satisfazem a Eq, (18), foram instáveis, conforme esperado.

Ao contrário do que foi observado no esquema avançado em v e z por Schleicher et al. (2003), a implementação direta não apresentou problemas. Pelo contrário, neste caso o esquema reverso fez com que os ruídos crescessem exponencialmente ocultando assim o pulso que continua propagando, ver Figura 2 (direita).

Também foram estudados numericamente os outros quatro esquemas de DF descritos neste trabalho. Não conseguimos implementar nenhum deles tal que produzisse resultados estáveis, nem em forma direta, nem reversa. Isto era de se esperar para o esquema avançado em v e centrado em z, que é incondicionalmente instável. Porém, a estabilidade prevista para os esquemas atrasado em v e avançado em z (para Δv > 0), atrasado em v e z (para Δv < 0) e atrasado em v e centrado em z (incondicionalmente estável) não pôde ser reproduzida numericamente. Isto se deve ao fato de que aplicamos uma técnica de análise de estabilidade válida para equações diferenciais parciais com coeficientes constantes para uma equação diferencial parcial com coeficientes não constantes, Lembramos que todos os esquemas tratados são implícitos, porém foram implementados de forma explícita pelo fato de usarmos condições de contorno nulas. Portanto, as conclusões obtidas para o caso da implementação de esquemas implícitos são as mesmas. As nossas implementações implícitas confirmam esta afirmação.

Mudança de variável-Caso β=ln (v/v_o):

O objetivo desta seção é fazer uma comparação entre as imagens obtidas quando fazemos a remigração do refletor utilizando o esquema de DF avançado em v e z (Schleicheret al., 2003) e as imagens obtidas quando fazemos mudanças de variáveis na equação da onda imagem.

Para tal comparação foi simulado numericamente um experimento sísmico como descrito na subseção anterior para os casos de mudanças de variáveis. Para cada um dos casos dados pela Eq. (22), usamos o respectivo valor para suas condições de estabilidade. Exemplificamos, então, para o caso β=ln (v/v_o). Para os fins de implementação, usamos uma velocidade normalizada v' =v/(1000 m/s) para o cálculo deβ=ln (v/v_o).

Para este caso foi utilizado o valor de Δβ= 0,002, onde Δβ é menor do que o valor de

que é equivalente à 0,002343, dado pela condição de estabilidade da Eq. (20).

Fazendo a remigração a partir de uma velocidade migrada incorreta, usando adiscretização representada pela Eq. (11), obtemos a Figura 3 que representa a propagação do refletor plano até este atingir sua profundidade correta com velocidade 3000 m/s.

Omitimos aqui as figuras correspondentes às outras mudanças de variáveis, pois estas são muito similares à já apresentada na Figura 3.

Portanto, as mudanças de variáveis não alteraram o resultado da remigração. Desta forma, a mudança de variável não ajudou em solucionar o problema do estiramento do pulso da onda imagem. Por outro lado, nos mostrou formas mais econômicas para a execução do algoritmo, uma vez que as três imagens da Figura 3 são práticamente idênticas àquelas obtidas em Schleicher et ai. (2003) para o caso avançado em v e z reverso, mas foram obtidas com um custo numérico diferente, idem para todos os outros casos. Como visto na Figura 1, a mudança que nos proporciona vantagens é a mudança β=ln (v/v_o), pois seu tempo computacional é o menor de todos os três tipos estudados.

Refletor não plano

Para mostrarmos a qualidade da equação da onda imagem como método de remigração em profundidade, apresentamos um modelo sísmico que contém um refletor plano inclinado e uma sinclinal, que aparecem freqüentemente em dados reais. Ao fazermos a migração deste modelo obtemos uma seção sísmica migrada como mostrada na Figura 5(a). Observamos um evento conhecido como "gravata" que é gerado pelo refletor sinclinal. Do estudo teórico deste tipo de refletor sabemos que esta gravata não se dissolveu porque a migração foi feita com uma velocidade incorreta muito abaixo do valor correto. Ao fazermos a remigração desta seção sísmica, este evento deve tornar-se uma sinclinal quando atingida a velocidade correta. Este efeito pode ser observado pela propagação na Figura 5, onde vemos que na velocidade v =6000 m/s a profundidade correta de ambos os refletores foi atingida e a gravata se dissolveu por inteira.

CONCLUSÃO

Neste trabalho, demos continuidade à análise numérica da equação da onda imagem para remigração em profundidade, que foi iniciado em Schleicher et ai. (2003). Estudamos seis esquemas de DF para esta equação. Analisamos sua consistência e estabilidade, assim como sua dispersão e dissipação. Destas análises concluímos que todos estes esquemas são consistentes com a equação da onda imagem. Pela análise de von Neumann, quatro são condicionalmente estáveis, um é incondicionalmente estável e um instável.

Vimos também que mudanças de variáveis na equação da onda imagem não solucionam nenhum dos problemas numéricos, pois os resultados teóricos enuméricos se mostraram muito próximos dos resultados encontrados sem tais mudanças (Schleicher et al., 2003, 2004). No entanto, para o caso β=ln (v/v_o), tivemos um ganho no ponto de vista computacional. Este ganho é de grande interesse quando trabalhamos com esquemas implícitos, pois as matrizes não precisam ser atualizadas acada iteração na velocidade.

Como no caso da variável v, a estabilidade e dispersão mostram uma forte dependência do número de ondas. Onde os números de ondas presentes nos dados a serem remigrados são mais altos na horizontal do que na vertical, a estabilidade pode ser de difícil obtenção. Quanto à dispersão e dissipação numérica, concluímos que é impossível, para todos os casos estudados, que estas sejam reduzidas para o problema inteiro, visto que as condições dependem da profundidade e da velocidade de migração. Portanto, estes efeitos podem somente ser reduzidos para níveis aceitáveis pela escolha de intervalos da malha muito pequenos.

Mediante um estudo numérico da equação da onda imagem, aqui representado por um modelo simples com um único refletor plano (outros tipos de refletores mais complexos são apresentados em Schleicher et ai. (2004)), vimos que imagens remigradas de refletores propagam à profundidade correta quando a velocidade atinge o seu valor correto e estruturas de gravatas se dissolvem corretamente. Estas aplicações numéricas indicam o potencial da remigração da onda imagem de ser útil como uma ferramenta para análise de velocidade de migração quando se possui informação adicional sobre a posição do refletor, por exemplo, a sua profundidade medida em um poço.

A estabilidade dos esquemas de DF foi testada nos experimentos numéricos. Para três dos esquemas estudados, os resultados numéricos confirmam as previsões teóricas. Para os outros três, não conseguimos produzir resultados estáveis conforme teoricamente previstos. Os esquemas avançados na variável da velocidade e avançados ou atrasados em profundidade foram os que tiveram os resultados mais proveitosos, pois estes esquemas satisfazem as previsões teóricas quanto à sua estabilidade e podem ser utilizados para realizar a propagação da onda imagem, sendo o primeiro para velocidades crescentes e o segundo para velocidades decrescentes.

AGRADECIMENTOS

Agradecemos aos dois revisores pelas sugestões que ajudaram a melhorar este trabalho, bem como à Petrobras, à FAPESP, ao CNPq e aos patrocinadores do consórcio Wave Inversion Technology (WIT) pelo suporte financeiro.

Recebido em 14 janeiro, 2008 / Aceito em 15 setembro, 2008

Received on January 14, 2008 / Accepted on September 15, 2008

NOTAS SOBRE OS AUTORES

Fernando Perin Munerato. Bacharel em Matemática Aplicada e Computacional pela UNICAMP em 2003. Mestre em Matemática Aplicada com ênfase em Geofísica Matemática pela UNICAMP em 2006. Estudante de doutorado em Ciências e Engenharia de Petróleo pela Faculdade de Engenharia Mecânica da UNICAMP. Seus interesses de pesquisa incluem métodos de imageamento sísmico, bem como simulação numérica de reservatórios de hidrocarbonetos.

Jõrg Schleicher. Graduado e mestre em Física pela Universidade de Karlsruhe - Alemanha , em 1985 e1000, respectivamente. Fez doutorado em Geofísica Aplicada no Instituto de Geofísica da Universidade de Karlsruhe. Em 1993, foi contratado pela mesma como colaborador científico. Desde 1995, trabalha na UNICAMP como professor do Depto, de Matemática Aplicada. Seus interesses de pesquisa incluem todos os métodos de modelagem, imageamento e inversão sísmicos. Em 1995, recebal o prêmio J. Clarence Karcher Award da Society of Exploration Gecphysicist.s

Amélia Novais. Bacharel em Matemática pela Universidade Federal Fluminense (UFF) em 1990, Mestre em Matemática pelo IMPA em 1993, Doutora em Matemática Aplicada com ênfase em Geofísica Matemática pelo IMECC/UNICAMP em 1998. Estagiou na Universidade de Karlsruhe, Alemanha , em 1997. Foi professora na Universidade Federal de São Carlos no período de 1996 a 2002, e desde então é professora do Departamento de Matemática Aplicada na UNICAMP. Seu interesse em pesquisa é em métodos sísmicos

• FOMEL S 1994. Method of velocity continuation ln the problem of Seilsmic time migration. Russian Geology and Geophysics, 35 100-111.
• FOMEL S 2003. Time migration velocity analysis by velocity continuation. Geophysics, 68: 1662-1672 .
• HUBRAL P, TYGEL M & SCHLEICHER J 1996. Seismic Image waves Geophysical Journal lnternational , 125 431-442.
• JAYA M, BOTELHO M, HUBRAL P& LlEBHARDT G. 1999. Remigration of ground-penetrating radar data. Journal of Applled Geophysics, 41 19-30.
• NOVAIS A, COSTA JC & SCHLEICHER J 2008. GPR velocity determination by image-wave propagation Journal of Applled Geophysics, 65 65-72.
• SCHLEICHER J, NOVAIS A & MUNERATO F. 2003. Análise numérica da equação da onda imagem na profundidade. TEMA -Tendências em Matemática Aplicada e Computacional, 4 267-276.
• SCHLEICHER J, NOVAIS A & MUNE RATO F. 2004. Migration velocity analysis by depth Image-wave remlgratlon Ilrst results. Geophyslcal Prospecting, 52 559-574.
• SCHLEICHER J, COSTA JC & NOVAIS A. 2008. Time migration velocity analysis by image-wave propagation of common-image gathers. Geophysics, 73 VE161-VE171
• STRIKWERDA J 1989. Finite Difference Schemes and Partial Differential Equations. Wadsworth & Brooks, California. 386 p.
• THOMAS J 1995. Numencal Partial Differential Equations Springer-Verlag, New York. 436 p.
• TYGEL M, SCHLEICHER J & HUBRAL P1994. Pulse distortion ln depth migration Geophysics, 59 1561-1569.

Datas de Publicação

• Publicação nesta coleção
  07 Maio 2009
• Data do Fascículo
  Dez 2008

Histórico

• Recebido
  14 Jan 2008
• Aceito
  15 Set 2008