Transformada de Laplace na solução de problemas inversos dinâmicos da sísmica

Transformada de Laplace na solução de problemas inversos dinâmicos da sísmica

Georgy Mitrofanov^I; Viatcheslav Ivanovich Priimenko^II; Roseane Marchezi Misságia^II; Luis Henrique Amaral^III

^IInstitute of Geology and Geophysics, Siberian Branch of the Russian Academy of Sciences, pr. Koptyuga, 4, Akademgorodok, 630090, Novosibirsk, Russia. Phone: [phone] - E-mail: [e-mail]

^IILaboratório de Engenharia e Exploração de Petróleo (LENEP), Universidade Estadual do Norte Fluminense Darcy Ribeiro (UENF), Rod. Amaral Peixoto, km 163,Imboassica, 27925-310 Macaé, RJ, Brasil. Tel.: [phone]; Fax: [phone] - E-mails: [e-mail]; [e-mail]

^IIIPetrobras/Exploração/Geofísica Aplicada à Exploração/Processamento Geofísico, Av. República do Chile, 65, Centro, 20031-912 Rio de Janeiro, RJ, Brasil. Tel.: [phone] - E-mail: [e-mail]

RESUMO

Neste artigo são analisados, utilizando métodos de otimização, vários aspectos da aplicação das transformadas de Laplace temporal e de Fourier-Bessel espacial na solução de problemas inversos dinâmicos da sísmica. Modelos delgados do meio são considerados, como exemplo, para demonstrar que a parte real do parâmetro de Laplace pode ser utilizada como um regularizador nas duas principais etapas de construção da solução do problema inverso: solução do problema direto e combinação da solução teórica do problema direto com os espectros dos traços sísmicos. Isto permite evidenciar que a partir da variação deste parâmetro é possível melhorar as propriedades do funcional de erro, acelerar a convergência do método e obter uma aproximação inicial mais acurada. Além disto, estudou-se a influência do parâmetro de Laplace e dos filtros de suavização no grau de similaridade entre os dois tipos de espectros: calculado usando os sismogramas (considerando a limitação da abertura de observações reais) e obtido com base nas fórmulas analíticas. Os resultados deste estudo possibilitaram o desenvolvimento de procedimentos computacionais, que visam garantir a boa qualidade do cálculo dos espectros discretos, usando para isto as transformadas de Fourier-Bessel e de Laplace.

Palavras-chave: sistema de Lamé, meio estratificado, reservatório delgado, domínio espectral, sistema de observação multi-onda e multicomponente, problema inverso dinâmico.

ABSTRACT

In this article we analyze several aspects of using the temporal Laplace and spatial Fourier-Bessel transforms in the solution of seismic inverse dynamic problems by optimization methods. Considering, as an example, thin layer models of media, we show that the real part of the Laplace parameter can be used as a regularization parameter in two principal stages of the solution of the inverse problem: solution of the direct problem, and combination of the theoretical solution of the direct problem with spectra of seismograms. It is shown that varying this parameter we can improve the properties of the misfit functional, accelerate the convergence of the method, and get a better initial approximation. Besides this, there was studied the influence of the parameter of Laplace and smoothing filters in the degree of similarity between two types of spectra: calculated using the seismograms (considering the limitation of the aperture of real observations) and obtained on basis of the analytical formulas. The results of this study made possible the development of computational proceedings, which aim to guarantee the good quality of the calculation of discrete spectra, using for this the Fourier-Bessel and Laplace transforms.

Keywords: Lamé system, stratified media, thin layer reservoir, spectral domain, multi-wave and multi component observation, inverse dynamic problem.

INTRODUÇÃO

A propagação de ondas sísmicas no meio rochoso é um processo complexo, caracterizado por vários parâmetros, entre eles podemos citar os cinemáticos e dinâmicos. Os parâmetros cinemáticos referem-se ao estudo de frente de ondas e raios, e estão estreitamente relacionados ao cálculo do tempo de percurso das ondas. Os parâmetros dinâmicos incluem amplitude e energia, forma da onda e características espectrais, polarização espacial, e peculiaridades de interferência. Esta divisão entre parâmetros cinemáticos e dinâmicos é convencional, especialmente, se considerarmos que na solução de problemas diretos e inversos, e, também, na supressão dos ruídos presentes no sinal, é necessário utilizar as características cinemáticas das ondas. Ao mesmo tempo, é razoável considerar a cinemática das ondas como um objeto independente da geofísica teórica, visto que importantes características do meio, tais como: profundidade, forma dos refletores sísmicos, velocidades de propagação das ondas, são medidas com maior precisão através da análise dos tempos de percurso das ondas, do que, por exemplo, através das amplitudes das ondas registradas.

A solução de problemas inversos para definição dos parâmetros do objeto alvo, usando o campo sísmico registrado, se configura como uma das etapas mais importantes do processamento e interpretação de dados sísmicos. Em função do objeto investigado podemos dividir tais problemas nas seguintes classes:

1) Determinação dos parâmetros da fonte do campo de onda sísmica, supondo que as características do meio sãoconhecidas.

2) Determinação das características do meio, supondo que os parâmetros da fonte são conhecidos.

Os problemas do primeiro tipo são usualmente encontrados na sismologia, na definição dos parâmetros de epicentros dos terremotos. Os problemas do segundo tipo são típicos na exploração sísmica, por exemplo, na determinação da velocidade e densidade de um modelo estratificado do meio.

A depender dos parâmetros do campo de onda usados na solução dos problemas inversos é possível distinguir entre problemas cinemáticos e dinâmicos. Nos problemas inversos cinemáticos o tempo de trânsito das ondas sísmicas é utilizado como uma informação adicional. A solução dos problemas inversos dinâmicos é caracterizada pela utilização da informação sobre a amplitude das ondas sísmicas.

Atualmente os métodos de solução dos problemas inversos cinemáticos ainda são mais desenvolvidos. Alekseev et al. (1979), Nolet et al. (1987), Cervený (2001) abordam as formulações gerais de tais problemas, condições de existência das soluções, questões de estabilidade e unicidade, e, também, propõem algoritmos para solução. Recentemente a teoria dos problemas inversos dinâmicos sísmicos vem se desenvolvendo intensivamente, veja Alekseev (1967), Romanov (2002), Lavrentiev et al. (2003). No entanto é importante ressaltar, que até para o problema 1D, a utilização da solução obtida é limitada em virtude da idealização dos modelos, e, também, devido a dificuldades computacionais. Os resultados obtidos são pouco expressivos, e justificados pelo baixo desenvolvimento da teoria dinâmica da propagação das ondas elásticas no meio em subsuperfície; pelo baixo conhecimento da distribuição dos parâmetros no meio geológico; pela dependência entre as características dinâmicas das ondas e as propriedades físicas das rochas. Desta forma, constata-se que apesar da história de desenvolvimento da teoria dos problemas dinâmicos inversos aplicados a sísmica, até o presente momento, os métodos não são suficientemente efetivos para serem utilizados com êxito na exploração sísmica. Os métodos de otimização da solução dos problemas inversos, baseados na solução aproximada dos problemas diretos, particularmente, o método de traçamento de raio, são amplamente ilustrados, por exemplo, em Kendall & Stuart (1973), Jurado et al. (1995). Porém, nestes trabalhos, o problema é considerado, sob o ponto de vista da sísmica tradicional, para meios geológicos com camadas espessas, onde o tempo de trânsito da onda até as camadas supera metade do comprimento da onda. No caso de aplicação desta concepção para modelos de camadas delgadas, varias questões ficaram sem resposta, tais como: o quanto o campo de onda é adequadamente descrito, quais tipos de ondas devem ser levados em conta, entre outras.

Muitos trabalhos teóricos foram dedicados à solução dos problemas inversos em modelos delgados e verticalmente heterogêneos, ver, por exemplo, Alekseev (1967), Blagoveshchenskii (2001), Kabanikhin & Lorenzi (1999), Romanov (2002). Estesautores construíram algoritmos teóricos e numéricos que comprovam a possibilidade de determinação dos coeficientes de Lamé e densidade em função da profundidade. Mas em tais formulações matemáticas foram usadas suposições sobre o registro do campo elástico em qualquer ponto da superfície livre, e determinados tipos de fontes, que permitem reduzir o sistema de equações elásticas a um sistema mais simples, transformando o problema inverso completo em uma série de problemas inversos, relativamente simples. Porém, tais fontes são pouco aplicadas na indústria, e sua criação requer gastos excessivos. Em virtude disto, investigamos o problema inverso 1D, utilizando uma fonte do tipo centro de dilatação instantânea, que serve como modelo de fonte explosiva na terra e na água. Neste caso, a redução do sistema de equações elásticas em equações escalares separadas torna-se inviável, requerendo a investigação da solução do problema inverso completo, o que representa laborioso problema matemático e numérico. Neste artigo a solução deste problema é alcançada através dos métodos de programação não-linear, minimizando o funcional objetivo, construído no domínio espectral. Com isso, construímos a solução exata do problema direto considerando um sistema de elasticidade relativamente complexo para camadas delgadas. Usamos para isso, uma concepção matricial de redução do sistema de elasticidade ao sistema de equações do tipo Riccatti,proposto por Thomson e Haskel (Thomson, 1950; Haskel, 1953), veja também Ursin (1983) e Molotkov (1984).

As investigações prévias do funcional de erro, executadas no domínio espectral, mostraram a sua não convexidade (Karchevsky, 2005) e existência de máximos e mínimos locais. Tais resultados representam uma conseqüência direta da não-linearidade do problema inverso, e presença de muitos parâmetros (freqüência espacial, freqüência temporal, quantidade de freqüências espacial e temporal ao construir funcional, e combinação destes parâmetros durante a solução numérica do problema inverso). Mas, mesmo assim, estas investigações, realizadas a partir de exemplos numéricos, proporcionaram a formulação de uma estratégia de minimização, que leva ao mínimo global. A implementação desta estratégia permitiu estimar os parâmetros do reservatório delgado, cujas espessuras das camadas são similares às dos reservatórios reais. Na construção de um funcional de erro no domínio espectral surge um problema com a apresentação do campo de onda espalhado, registrado na superfície, usando transformadas integrais. No caso do modelo axialmente simétrico, esta representação faz uso das transformadas de Laplace com respeito ao tempo, e de Fourier-Bessel com respeito às coordenadas espaciais. Note a importância do estudo das características especiais a partir do uso das transformações indicadas no cálculo dos espetros bidimensionais, e na construção da solução de problemas inversos no domínio espectral, quando se utilizam as ondas espalhadas, registradas na superfície livre.

Neste artigo apresenta-se a investigação das principais possibilidades do método proposto.

FORMULAÇÃO DE PROBLEMA

Problema direto

Consideremos um meio composto por n camadas estratificadas: 0 = z₀ < z₁ < ... < z_n < ∞. As propriedades físicas de cada camada z_k-1 < z < z_k, k = 1,2, ..., n, caracterizam-se pelos coeficientes de Lamé λ,m, e densidade r - funções constantes por partes com descontinuidades nos pontos z_k, k = 1,2,..., n. As oscilações elásticas no meio são geradas por uma fonte do tipo centro de expansão, definidopela fórmula:

ρ_sπgrad δ(x, y, z-z_s)g(t)

onde z_s (z_s > 0, z_s≠ z_k, k = 1,...,n) é a profundidade da fonte, g(t), g(t) º 0 para t < 0, é a forma de impulso, δ(x,y,z) é a função de Dirac, e ρ_s - a densidade da camada onde a fonte está situada.

Visto que o problema é formulado para um meio horizontalmente estratificado, podemos discorrer sobre a simetria do meio e fonte, e representar o problema no sistema de coordenadas cilíndricas:

x = r cos θ, y = r sin θ, z = z

0 < r < ∞ 0 < θ < 2π, - ∞ < z < ∞

Neste caso o vetor de deslocamento não depende da coordenada angular e tem somente duas componentes: u(z,r,t) - radial (horizontal) e ω(z,r,t) - vertical. O sistema de Lamé é representado em coordenadas cilíndricas, da seguinte forma:

A transição de uma função s(z,r,t) para o domínio espectral é realizada utilizando as transformadas de Fourier-Bessel com respeito à variável espacial r, e Laplace com respeito à variável temporal t:

onde J_m é a função de Bessel da ordem m, ν é a freqüênciaespacial da transformada de Fourier-Bessel, e p é o parâmetro da transformada de Laplace; p = -α+ιω e ω = 2πƒ, ƒ é a freqüência temporal. O índice m é igual a 0 para a componente vertical w(z,r,t), e 1 para a componente horizontal u(z,r,t).

As Eqs. (1) podem ser representadas no domínio espectral como:

onde (p) é a transformada de Laplace da função g(t). Como o problema direto, considera-se a definição das funções (z,n,p) (z,n,p), que no domínio z > 0 satisfazem as Eqs. (3) e as seguintes condições de contornos e de atenuação:

Além disto, nos pontos de descontinuidade dos coeficientes são válidas as seguintes condições de equilibro:

= ƒ(z+0) - ƒ(z-0) representa o salto da função ƒ no ponto z.

As condições definidas pelas Eqs. (3)-(5) são típicas da formulação de problemas matemáticos aplicados à exploração sísmica.

Problema inverso

O problema inverso consiste na definição das velocidades V_p = , densidade r, e espessura de cada camada, usando como informação adicional a solução do problema direto, definido pelas Eqs. (3)-(5):

onde

^R (ν,p), ^R (ν, p), ν ∈ Ω_v, ω ∈ Ω_ω são funções conhecidas,

e N_ν, N_ω são números finitos.

Numericamente o problema inverso, definido pelas Eqs. (3)-(6), pode ser resolvido minimizando o seguinte funcional objetivo:

que representa uma função dos parâmetros objetivos Θ (características elásticas V_p, V_s, ρ e espessura das camadas). Na construção deste funcional é preciso transformar os sismogramas em espectros bidimensionais ^R (ν,p), ^R(ν,p), e também calcular os valores teóricos destes espectros, utilizando para isto, a solução do problema direto no domínio espectral (0, ν, p, Θ), (0, ν, p,Θ). As constantes h_ν, h_ω são os multiplicadores de normalização, e dependem da quantidade das freqüências espacial e temporal utilizadas, e do intervalode amostragem:

Na Eq. (7), os índices k, n foram omitidos para simplificação.

INFLUÊNCIA DAS TRANSFORMADAS NO FUNCIONAL J(Θ) E ESPECTROS CALCULADOS

Propriedades do funcional objetivo

O funcional J(Θ) depende tanto do vetor Θ como dos parâmetros h_ν, h_ω, Ω_ν, Ω_ω, a. A solução do problema de minimização é reduzida a uma busca de todos, ou uma parte dos componentes do vetor Θ, que minimizem a Eq. (7). Note que a minimização do funcional objetivo no problema inverso é um dos métodos mais divulgados, ver, por exemplo, Himmelblau (1972), Gill et al. (1981). A dificuldade principal de minimização deste funcional reside na existência de mínimos e máximos locais. Posteriormente será mostrado que o funcional objetivo, definido pela Eq. (7), também apresenta esta dificuldade.

Vamos analisar as mudanças do funcional objetivo em relação aos diferentes parâmetros. Visto que os valores de J(Θ) são calculados utilizando as funções (z,ν,p), (z,ν,p) é importante compreender o comportamento destas funções com respeito a alguns dos principais parâmetros. Os experimentos, realizados no cálculo dos espectros bidimensionais para vários modelos, mostraram que os parâmetros a, n têm maior influência na estrutura destas funções e do funcional, também.

A investigação da solução do problema direto possibilitou a obtenção das características espectrais das funções , .A principal conclusão de tais investigações é que para α > 5 os valores das funções , praticamente coincidem para modelos com escalas variáveis de freqüências espacial e temporal. A partir da redução gradual de α, as diferenças nos valores das funções , começam a se manifestar para os meios com camadas relativamente espessas. Para os valores mais baixos deste parâmetro, α < 0.1, as diferenças tornam-se essenciais para todos os modelos delgados utilizados. O sentido físico deste efeito pode ser entendido a partir da análise da transformada de Laplace. Esta transformada pode ser interpretada como a transformada de Fourier do sismograma original com função de peso e^-αt. Por isso, para valores elevados de α, a contribuição de tempos, correspondentes a reflexões de horizontes mais profundos e outros tipos de ondas (múltiplas, convertidas, entre outras), nos valores das funções , serão praticamente nula. Estas reflexões serão anuladas pelo fator exponencial. E ao contrário,reduzindo os valores de α, aumentamos a contribuição destes componentes nos valores das funções , .

A partir da análise da estrutura da Eq. (7), vários experimentos numéricos, alguns bastante complicados, foram executados para uma sistematização. Neste artigo, apresentaremos somente um exemplo para descrever as regras definidas durante a investigação. Os resultados deste experimento são ilustrados na Figura 1, e apresenta a estrutura do funcional J(Θ) para o modelo definido na Figura 2, considerando a variação do parâmetro V_p na última camada do modelo - camada 4. A camada e o parâmetro foram escolhidos de forma arbitrária. Durante a construção do funcional foram usados somente os valores da função (0, ν, p) para a representação do contraste de suas singularidades. É importante pontuar que as regularidades básicas encontradas no comportamento de funcional foram similares para outros modelos delgados e combinações das funções , .

É óbvio que variando mais parâmetros em várias camadas, o comportamento do funcional objetivo torna-se mais complicado. Porém, este simples exemplo, possibilita a análise de algumas singularidades e formular as recomendações para elaboração de uma estratégia de minimização deste funcional. Na Figura 1 o eixo horizontal corresponde a , e o eixo vertical - aos valores do funcional J(Θ). Esta investigação tem como meta analisar a influência das freqüências temporal e espacial no comportamento do funcional construído.

A Figura 1(a) mostra o comportamento do funcional para as freqüências localizadas numa área, situada na parte esquerda superior da figura. Neste experimento, utilizamos o seguinte conjunto de parâmetros: α = 0.1, ν = {0,0.1}, e 100 valores de ƒ regularmente distribuídas no intervalo [5, 50] Hz.Note que o comportamento do funcional, construído nesta área, não é favorável à busca do mínimo global, observa-se um grande número de extremos locais irregulares. Isto impossibilitou a definição do valor exato = 4.84.

A Figura 1(b) mostra o comportamento do funcional com relação ao aumento do intervalo e quantidade de freqüências temporais. Neste caso α = 3 e ν = {0,0.01}. Os gráficos apresentados nesta figura correspondem aos seguintes intervalos de freqüências temporais: [5, 50], [5, 100], [5, 150], [5, 200], [5, 300], [5, 500], [5, 900], [5, 1900], todos com intervalo de amostragem de 1 Hz. O primeiro gráfico inferior (mais inclinado) corresponde ao intervalo mínimo de freqüências, o último gráfico superior (mais abrupto) corresponde ao intervalo máximo de freqüências. É visível a tendência de suavização do funcional e desaparecimento dos mínimos locais. Podemos observar que a partir de alguns intervalos este efeito começa a desaparecer. A Figura 1(c) ilustra o comportamento do funcional em resposta ao aumento do número das freqüências temporais em um intervalo limitado. Aqui α = 3, ν = {0, 0.01}, o intervalo das freqüências [5, 50] Hz, e os gráficos do funcional J(Θ) foram construídos para as seguintes partições deste intervalo:10, 50, 100, 500, 2000. O gráfico mais inferior corresponde à partição mínima (os três últimos gráficos coincidiram entre si, praticamente). Observe a suavização dos gráficos para velocidades menores, quando < 2. A Figura 1(d) mostra a variação do funcional em relação ao aumento do número das freqüências temporais em um intervalo limitado. Aqui foram utilizados os seguintes parâmetros: α = 3, 100 valores da ƒ no intervalo [5, 50] Hz. Para cada um dos três gráficos, 2, 11, e 21 valores de ν amostrados com regularidade foram escolhidos no intervalo [0, 0.01]. Observe que o aumento do número das freqüências espaciais praticamente não muda o comportamento do funcional.

É interessante analisar o comportamento do funcional em relação às baixas freqüências temporais, veja a Figura 1(e). Onde α = 3, ν = {0,0.01}. Para a construção dos cinco gráficos foram escolhidos os seguintes intervalos das freqüências temporais: [5, 10], [5, 20], [5, 30], [5, 40], [5, 50]. Para cada um dos intervalos foram escolhidos 30 valores de freqüências, amostrados com regularidade. O gráfico mais inclinado corresponde ao primeiro intervalo, e o mais abrupto corresponde ao último. Isto possibilita observar com clareza que o funcional objetivo tem uma forma quase convexa nos intervalos de baixas freqüências, quando ƒ < 20 Hz. Nas freqüências médias, o incremento deste intervalo acentua o declive de curva e aparecem os extremos locais. A Figura 1(f) ilustra o comportamento do funcional em relação ao parâmetro α. Onde ν = {0,0.01}, e foram utilizados 100 valores da ƒ no intervalo [5, 50] Hz. Na construção destes gráficos usamos os valores de α iguais a 1 (gráfico inferior) e 3 (gráfico superior) - compare com o gráfico da Figura 1(a), quando α = 0.1. Note que o crescimento de α provoca uma mudança no comportamento do funcional para uma forma mais convexa e suave.

Os resultados obtidos permitem tecer conclusões que são úteis na solução numérica do problema inverso, baseada na minimização do funcional objetivo.

1. O crescimento do parâmetro a suaviza os espectros bidimensionais (z, ν, p), (z, ν, p).

2. Para elevados valores de a meios com camadas de grande espessura e meios com estrutura mais complexa são praticamente indistinguíveis.

3. Com ampliação do intervalo das freqüências temporais utilizadas e aumento do intervalo de amostragem o funcional objetivo mostra uma tendência à suavização, Figuras 1(b) e 1(c).

4. A busca do ponto de mínimo global do funcional torna-se impossível se utilizarmos os parâmetros, para quais os valores das funções (z, ν, p), (z, ν, p) pertençam ao domínio onde a função de resposta do meio apresenta "fortes" oscilações, relacionadas ao espectro do problema direto, Figura 1(a).

5. Aumento do número das freqüências espaciais utilizadas em um intervalo limitado não melhora as propriedades do funcional objetivo, Figura 1(d).

6. O funcional objetivo tem uma forma quase convexa e mais suave para as baixas freqüências temporais. A ampliação do intervalo das freqüências temporais, com inclusão das freqüências médias e altas, transforma o funcional num mais abrupto. Mas neste caso, aparecem os extremos locais, conforme ilustra a Figura 1(e).

7. Com aumento de a o funcional objetivo transforma-se em um funcional quase convexo e mais suave, e ao contrário, com a diminuição do parâmetro α o funcional transforma-se num mais abrupto. Mas neste caso, aparecem os extremos locais, conforme mostra a Figura 1(f).

A utilização de tais conclusões melhora a seleção dos parâmetros usados na construção de J(Θ); em particular, dos parâmetros α e ν, dos intervalos das freqüências temporal e espacial, e do número destas freqüências. A melhor escolha de tais parâmetros possibilita formar uma estrutura do funcional objetivo favorável à minimização. A mudança deste conjunto, em algumas etapas da minimização, permite melhorar a busca do mínimo global e, por esta razão, garantir a solução do problema inversocom maior precisão.

ESTRATÉGIA DE MINIMIZAÇÃO DO FUNCIONAL J(Θ) USANDO O PARÂMETRO α

Como mostrado acima, o funcional objetivo pode ter uma estrutura complexa, incluindo extremos locais. É conhecido, que a busca do mínimo global é um problema extremamente difícil de ser resolvido. Portanto, a escolha de uma estratégia de minimização de tal funcional é uma etapa importante na solução geral. A melhoria das propriedades do funcional e simplificação do processo de minimização poderia ser feita estendendo o intervalo das freqüências temporais. Porém, a presença no sinal sísmico real de uma faixa de freqüências muito limitada, praticamente inviabiliza este caminho. Por isto, o mais prático é um esquema baseado nos diferentes valores do parâmetro α.

A estrutura deste esquema é a seguinte:

• Começar a minimização do funcional escolhendo como aproximação inicial um modelo composto de algumas camadas espessas;

• Nas primeiras etapas da minimização do funcional usar os maiores valores de α e intervalo de baixas freqüências temporais. Quando a velocidade de minimização do funcional começar a diminuir, podemos ampliar o intervalo de freqüências temporais e continuar o procedimento de minimização;

• Após a redução de α, mais uma vez implementa-se a minimização. Desta forma, reduzimos α por etapas, até quando o valor do parâmetro buscado da próxima camada espessa permanecer inalterado;

• Representamos a próxima camada espessa como umpacote de camadas mais finas. Inicialmente, calculamos a aproximação seguinte em freqüências baixas, depois disto, indicamos esta aproximação com exatidão, utilizando um intervalo de freqüências temporais mais amplo;

• Deste modo, passamos gradualmente a reconstrução dos parâmetros das camadas superiores para camadas inferiores.

Esta estratégia permite, utilizando um intervalo com as freqüências baixas e um valor elevado do parâmetro α, chegar mais perto do mínimo global. Para este conjunto de parâmetros o funcional tem uma forma em declive. Também é conhecido que para funcionais deste tipo, os métodos de gradientes apresentam um "baixo" desempenho na vizinhança do ponto mínimo, veja Nocedal & Wright (2006). Mas, esta aproximação chega mais perto do mínimo global que a aproximação inicial. Em seguida, passo a passo, diminuindo α e estendendo o intervalo de freqüências temporais, chegamos mais perto do mínimo global. A Figura 1(e) mostra que os extremos locais surgem, inicialmente, nos domínios de freqüências muito altas e muito baixas, e depois, com a ampliação do intervalo de freqüências temporais aproximamos mais do mínimo global. Se diminuirmos α no funcional objetivo pode ser observado um comportamento similar. Assim, selecionando o intervalo de freqüências temporais e variando a, é possível aproximar-se do mínimo global, e evitar os extremos locais. Na realização prática deste esquema, para cada pacote de camadas, foram utilizados os seguintes passos, que definem o algoritmo da solução do problema inverso no domínio espectral:

1. Reconstrução da velocidade das ondas longitudinais V_p:

• Escolha de

  ν = 0;

• Escolha de uma aproximação inicial

  V_p = const, em todas as camadas do pacote dado;

• Escolha do intervalo de freqüências [10, 25] Hz;

• Minimização do funcional objetivo, parando o processo quando funcional ou norma do gradiente do funcional for menor que um número positivo ε

  ₁;

• Escolha do intervalo de freqüências [10, 85] Hz;

• Utilização da aproximação obtida

  V_p, como uma aproximação inicial, para todas as camadas deste pacote, e minimização do funcional objetivo para as camadas mais finas;

• Processo de minimização do funcional deve ser interrompido quando o funcional ou a norma do gradiente do funcional for menor que um número positivo ε

  ₂.

2. Reconstrução da velocidade das ondas transversais V_s:

• Escolha de

  ν no intervalo [0.001, 0.01];

• Escolha de uma aproximação inicial

  V_s = c₁V_p+c₂ (

  c₁,c₂ são as constantes correlativas, típicas para área dada), utilizando os valores reconstruídos para todas as camadas do pacote analisado;

• Escolha do intervalo de freqüências [5, 65] Hz;

• Minimização do funcional objetivo, interrompendo o processo quando o funcional ou a norma do gradiente do funcional tornar-se menor que um número positivo ε

  ₃.

3. Escolha de um ou vários valores para ν do intervalo [0.001, 0.01], e minimização do funcional objetivo para a precisão de V_p, V_s.

A Figura 3 representa um exemplo da realização deste esquema para a definição do parâmetro V_p. Neste caso, inicialmente considerou-se o meio como homogêneo, contendo uma camada com 100 m de espessura na parte superior do modelo, veja a Figura 3(a). A fonte pontual foi posicionada na superfície livre z = 0. A seguir definem-se os parâmetros intervalares da primeira camada espessa, veja a Figura 3(b). Nesta etapa foram utilizados α = 10 e um intervalo de freqüências temporais relativamente baixos de [5, 15] Hz. Por força das condições impostas aos parâmetros do meio, os parâmetros da camada superior foram definidos com exatidão. Depois disto, de acordo com o esquema geral, definimos os parâmetros do modelo delgado dentro desta camada. Utilizaram-se α = 5 e o intervalo de freqüências temporais de [5, 25] Hz. Os resultados obtidos estão representados na Figura 3(c). Em todas as camadas deste conjunto os valores estimados correspondem aos verdadeiros, exceto ascamadas 9 e 10, onde visualizamos alguns erros.

A seguir passamos à próxima camada espessa, veja Figura 3(d), correspondente à segunda camada do modelo. Os parâmetros desta camada foram estimados e os parâmetros do pacote anterior, de camadas finas, foram indicados com exatidão. Depois disto, estimamos os valores das velocidades do pacote de camadas finas, correspondente à segunda camada espessa,Figura 3(e). Para isto, utilizamos α = 3 e intervalo de freqüências temporais de [5, 28] Hz, com partição em 47 freqüências. Isto resultou em uma boa aproximação para o modelo original. Os parâmetros são praticamente iguais aos valores verdadeiros, além das camadas 18, 19 e 20. Em seguida, passamos à terceira camada espessa e repetimos os passos executados nas camadas anteriores. Os resultados obtidos são apresentados nas Figuras 3(f)-3(h). Observe que todas as velocidades do modelo foram reconstruídas com alta precisão.

ALGUMAS PROPRIEDADES DAS TRANSFORMADASDE LAPLACE E DE FOURIER-BESSEL DISCRETAS

Para a solução do problema inverso é preciso representar os sismogramas no domínio espectral bidimensional usando as transformadas de Laplace e de Fourier-Bessel. Apesar desta operação não apresentar nenhuma dificuldade, é muito importante que sejam analisadas suas propriedades relacionadas à forma discreta dos dados, e a questão de limitação da abertura, com respeito às variáveis temporal e espacial. É necessário dizer que existem poucas publicações sobre esta questão. Diante disso, somente podemos tecer algumas conclusões sobre as propriedades discretas da transformada de Laplace, baseadas nas propriedades da transformada discreta de Fourier, visto que em relação à transformada de Fourier-Bessel esta questão praticamente não foi estudada.A investigação das propriedades discretas destas transformações será feita em base de funções analíticas. Neste caso a partir da análise dos resultados obtidos será possível examinar algumas propriedades elementares de análogos discretos das transformadas de Laplace e de Fourier-Bessel.

Consideremos uma função de variáveis x,t

onde θ(.) é uma função de Heaviside, definida pela fórmula

Assim, a função t₀(x) caracteriza o início do sinal e os parâmetros (α₀, α₁, ω₀) definem sua estrutura no sismograma, respectivamente. A escolha desta função está relacionada a algumas de suas propriedades. Primeiro, esta função permite, de uma forma muito elementar, considerar uma aproximação das propriedades das transformadas à estrutura dos sinais sísmicos observados. Segundo, para esta função podemos construir representações analíticas das transformadas de Laplace e de Fourier-Bessel. Terceiro, apesar de sua forma elementar é possível, usando estas representações analíticas, verificar a influência da intensidade inicial do sinal, ou alterações no tempo de trânsito do sinal resultante.

Consideremos dois casos da função t₀(x):

• I:

  t

  ₀(

  x) =

  t

  ₁ - const

• II:

  t

  ₀(

  x) =

  t

  ₁

  x

Estes casos possibilitam o entendimento de algumas propriedades das transformações discretas com respeito à mudança no tempo de trânsito do sinal.

Para representar a função U(x, t) no domínio espectral utilizamos as transformações de Laplace e Fourier-Bessel, definidas pela Eq. (2):

Relembramos que ν é a freqüência espacial da transformada de Fourier-Bessel, e p é o parâmetro de transformada de Laplace; p = -α+ iω e ω = 2πƒ, onde ƒ - é a freqüência temporal. O índice m é igual a 0 para a componente vertical do campo de onda, e 1 para a componente horizontal. Consideremos o caso quando o índice m é igual a 0 - um análogo da componente vertical de campo de onda.

Através da forma explicita da função U(x, t) no caso I obtemos a seguinte expressão do espectro bidimensional:

No caso II o espectro bidimensional completo da função U(x, t) será o seguinte:

onde , φ = arg(+ ν² - - 2iω₁a₁), a₁ = α₁ + αt₁, e ω₁ = ωt₁. Note, que para t₁ = 0 a Eq. (10) é equivalente à Eq. (11).

A imagem de Laplace da função U(x, t) pode ser representada na seguinte forma:

Agora analisaremos alguns resultados. Primeiro, consideremos a transformada de Laplace para o caso I, quando t₁ = 0 ou t₁ > 0. A estrutura e os parâmetros dos sinais gerados pela função são dados na Figura 4. Por conveniência, a função com α₀ = 100 representa um impulso curto, sendo definida como uma função do tipo I, e a função com α₀ = 10 representa um impulso longo, definido como uma função do tipo II. Os dados sintéticos, calculados através da Eq. (8), podem ser representados na forma de um sismograma. A geometria de aquisição fonte-receptor consiste de 01 família de tiro comum com 24 receptores espaçados de 50 m, sendo a distância entre a fonte e o primeiro receptor de 100 m. Assim define-se a abertura espacial de [100, 1250] m e temporal de [0, 3] seg. As Figuras 5 e 6 apresentam os resultados de alguns testes relacionados à transformada de Laplace. As partes real e imaginária da função foram obtidas a partir da Eq. (12), e estão representadasnas Figuras 5(a) e 6(a).

Os primeiros resultados, obtidos utilizando os procedimentos desenvolvidos, mostraram que já no nível da transformada de Laplace podem aparecer algumas singularidades computacionais, que deformam consideravelmente a estrutura dos espectros calculados. Por exemplo, durante o cálculo das partes real e imaginária de , caso t₁ = 0, surge uma deformação considerável da estrutura dos espectros calculados, principalmente no domínio de altas freqüências - comparar as Figuras 5(a) e 5(b).A Figura 5(a) representa as partes real e imaginária desta função, calculadas analiticamente, e a Figura 5(b) - calculadas numericamente. Através da análise dos intervalos mais restritos, é possível identificar que esta deformação é insignificante para freqüências menores que 50 Hz. Entretanto, a solução do problema inverso não pode ser restrita a freqüências menores que 50 Hz. Mas é óbvio, que para a solução do problema pode ser necessário considerar as freqüências maiores que 50 Hz. Por esta razão é preciso eliminar as componentes que provocam as deformações da estrutura dos espectros calculados. A nossa prática de utilização da transformada discreta de Fourier mostra que estas componentes, com tendência a crescer consideravelmente em altas freqüências, estão relacionadas aos efeitos laterais de "corte" dos impulsos do traço sísmico. A eliminação destes efeitos pode ser feita através da utilização de filtros especiais, ver Mitrofanov (1980). A utilização destes filtros, no caso da transformada de Laplace discreta, também garante uma melhoria significativa da qualidade dos espectros construídos, conforme mostra a Figura 5(c), onde são representadas as partes real e imaginária da transformada de Laplace da função U(x, t) do tipo I.

É possível notar que para a estrutura da função dada, caso t₁ > 0 (exemplo representado na Figura 4(b)), não é necessário utilizar os filtros especiais. Isto possibilita calcular os espectros de com base na transformada de Laplace discreta, que corresponde exatamente aos valores da , calculados analiticamente, conforme mostra a Figura 6. Nestas figuras, escolheu-se a faixa de freqüências até 60 Hz, para que fosse possível estimar com exatidão a semelhança entre os espectros construídos analiticamente e numericamente. Com relação aos resultados obtidos podemos ressaltar que a extensão temporal do impulso para t₁ = 0 provocou a compreensão dos resultados da transformada de Laplace e a aparência de um componente harmônico nas partes real e imaginária da função para t₁ > 0. Estesresultados coincidem com as propriedades conhecidas da transformada de Fourier.

Consideremos, agora, a análise dos espectros bidimensionais começando do caso t₁ = 0. A análise dos resultados obtidos será efetuada somente na componente real do espectro (a componente imaginária apresenta resultados semelhantes).A Figura 7(a) apresenta a parte real das componentes, correspondentes aos espectros bidimensionais dos sinais dos tipos I e II. Para o cálculo analítico dos espectros utilizamos a Eq. (10). Nestas figuras o eixo horizontal corresponde à freqüência espacial ν, e o eixo vertical à freqüência temporal ƒ.

Nas figuras apresentadas é mostrado o espectro correspondente a mudanças de ƒ no intervalo de 0.01 Hz até 60 Hz,para os dois tipos de sinais - I e II. Os resultados mostram que o cálculo dos espectros bidimensionais dos sismogramas, utilizando a integral dupla definida pela Eq. (9), provoca fortes distorções nas componentes espectrais. Além da tendência notada anteriormente nas freqüências altas, observa-se certa oscilação do espectro durante a variação da freqüência espacial, ver Figura 7(b). Ambos os efeitos podem ser retirados utilizando filtros. Para eliminação desta tendência no domínio das freqüências altas utilizamos filtros semelhantes aos utilizados no caso da transformada de Laplace, ver Figura 7(c). Para eliminar as variações adicionais falsas foram usados filtros mais suaves com respeito à variável espacial, veja a Figura 7(d). Estes filtros foram aplicados à transformada de Laplace discreta da função . A Figura 7(e) mostra a utilização dos filtros, seguida do procedimento desenvolvido neste artigo. Estes permitiram a obtenção de espectros bidimensionais, cujos valores coincidiram com os valores estimados de forma analítica a partir da Eq. (10), para ambos os tipos de sinais - I e II. Para t₁ > 0 é possível conseguir o mesmo efeito aplicando os filtros de suavização nos resultados da transformada de Laplace discreta. A não utilização destes filtros induz ao surgimento das mesmas oscilações no espectro durante a variação da freqüência espacial, notado no caso anterior, Figura 8(b). A aplicação destes filtros possibilita a obtenção de resultados semelhantes aos resultados obtidos a partir das fórmulas analíticos, Figura 8(c).

Todos os resultados apresentados demonstram a eficiência dos procedimentos desenvolvidos para a obtenção de sismogramas no domínio espectral. Por esta razão, alguns testes foram realizados com o objetivo de esclarecer melhor determinadas propriedades de tal transição. Um dos primeiros resultados pode ser a análise da influência do tempo de trânsito do sinal na estrutura dos espectros bidimensionais. Neste caso, isso é definido pela função t₀(x) e corresponde ao tipo II da função U(x, t).Note que as primeiras tentativas de encontrar correspondênciaentre o espectro analítico e o espectro construído a partir do desenvolvimento proposto nesta pesquisa, não produziram resultados positivos. Entretanto, uma investigação posterior revelou o motivo principal desta divergência entre os espectros: a integração com respeito à variável espacial x sobre um intervalo infinito permite "sentir" na estrutura do espectro bidimensional, definido pela Eq. (11), até as pequenas mudanças na função t₀(x). Contudo, se considerarmos o caso de abertura limitada é impossível identificar as pequenas variações. E, no caso da função do tipo II encontramos estes problemas.

As Figuras 9 e 10 representam os espectros bidimensionais, construídos utilizando os procedimentos desenvolvidos nesta pesquisa (Figs. 9(a) e 10(a)) e as fórmulas analíticas (Figs. 9(b) e 10(b)). Observe que para valores muito pequenos, t₁ < 0.00004 seg, ambos os espectros são praticamente iguais, conforme mostra a coluna esquerda das figuras. Porém, mesmo uma pequena variação desta grandeza, na ordem de 0.0001, é suficiente para provocar significativas diferenças entre os espectros. Para t₁ = 0.004 seg, percebe-se a sensibilidade do espectro, calculado a partir do sismograma, em relação ao tempo de trânsito do sinal para a variável espacial. Entretanto, o espectro construído a partir das fórmulas analíticas manifesta-se de uma forma bem diferente, conforme é ilustrado na coluna a direita das Figuras 9 e 10.

Para a solução do problema inverso é importante analisar a influência do parâmetro a na estrutura dos espectros obtidos. Uma das principais idéias da nossa pesquisa é utilizar o parâmetro a como um parâmetro regularizador para aproximar a solução exata do problema, obtida no domínio espectral, aos valores dos espectros bidimensionais, calculados utilizando os sismogramas observados. As Figuras 11 e 12 apresentam alguns resultados destes experimentos. A Figura 11 ilustra a influência do parâmetro α na estrutura da parte imaginária da função , calculada utilizando as fórmulas analíticas (coluna esquerda), e os procedimentos desenvolvidos nesta pesquisa (coluna direita). Os resultados apresentados nesta figura foram obtidos utilizando a função do tipo II com α₀ = 10 e t₁ = 0, valores do parâmetro α foram: 0.01(a), 10(b), 100(c). É evidente, que o aumento de α transforma a parte imaginária da função do tipo II numa função do tipo I (com α₀ = 100 e t₁ = 0). Este efeito foi observado através de uma comparação criteriosa entre os espectros correspondentes. O resultado fica evidente se levarmos em consideração como o parâmetro α é utilizado na Eq. (12). De acordo com esta fórmula o incremento de a é equivalente ao incremento do parâmetro do sinal α₀. Note que neste caso, no cálculo dos valores da transformada do Laplace discreta foram utilizadas todas as recomendações obtidas acima. Isto garantiu a similaridade entre os resultados obtidos analiticamentee numericamente.

A Figura 12 apresenta os resultados relacionados ao espectro completo, obtido com a utilização de ambas as transformações: de Fourier-Bessel e de Laplace. De acordo com estes resultados, o parâmetro α tem o mesmo efeito na variação do sinal, veja as Figuras 9 e 10 (com t₁ = 0.0004 seg). Isto pode ser justificado pela estrutura da Eq. (11). No entanto, esta influência será mais complexa, porque existe uma dependência entre α, α₁, t₁. Deste modo, a pode regularizar a estrutura do espectro bidimensional. Mas esta mudança influência na estrutura do impulso e não ajusta totalmente os espectros construídos a partir de fórmulas analíticas e os calculados com base nos análogos discretos das transformadas apontadas.

CONCLUSÕES

Na primeira fase do trabalho foi mostrado que através dos métodos de otimização é possível evidenciar que a parte real do parâmetro de Laplace permite melhorar essencialmente o processo da solução dos problemas inversos dinâmicos da sísmica. Este parâmetro controla a estrutura do funcional objetivo, transformando-o num funcional mais suave (ampliando essencialmente a escolha da aproximação inicial), ou num mais abrupto (acelerando a convergência do método utilizado).

Isto possibilita realizar uma estratégia eficaz de busca do mínimo global e da solução exata do problema. Entretanto, as vantagens indicadas relacionam-se somente à construção das soluções do problema inverso no domínio espectral. Mas os dados reais são observados no domínio espacial-temporal com uma substancial limitação do tempo de registro e da abertura. Assim, surge um grave problema de similaridade entre as observações reais e a solução teórica, construída no domínio espectral. Neste caso é necessário resolver diversas tarefas relacionadas com:

• Escolha de um modelo de referência (

  background model );

• Adaptação dos dados reais às condições do modelo teórico;

• Consideração nas soluções teóricas das particularidades dos dados reais, incluindo descontinuidade de dados sísmicos e limitação de abertura.

Na segunda parte do trabalho consideramos algumas das questões indicadas, usando as funções analíticas. Isto tornou possível conduzir os estudos das características especiais computacionais destas transformadas, e estudar a influência do parâmetro de Laplace e os filtros de suavização no grau de similaridade entre os dois tipos de espectros: calculado usando os sismogramas e obtido com base nas fórmulas analíticas. Os resultados deste estudo possibilitaram o desenvolvimento de procedimentos computacionais, que visam garantir a boa qualidade do cálculo dos espetros discretos bidimensionais, usando para isto as transformadas de Fourier-Bessel e de Laplace. Espera-se que este procedimento possibilite a criação de uma base para uma adaptação deste método, proposto para solução de problemas inversos no domínio espectral, para o caso dos espectros desinais reais.

AGRADECIMENTOS

Os autores agradecem aos revisores pelas valiosas sugestões e comentários, e à Petrobras pelo apoio dado para a execução deste trabalho.

Recebido em 18 Junho, 2008 / Aceito em 22 January, 2010

Received on June 18, 2008 / Accepted on January 22, 2010

NOTAS SOBRE OS AUTORES

Georgy Mitrofanov é graduado em Geologia e Geofísica pela Universidade Federal de Novosibirsk (UFN), Akademgorodok, Novosibirsk, Rússia, em 1972. Obteve seu doutorado em Física-Matemática em 1984 na UFN. Obteve seu título de Livre Docente em 1989 na UFN. É chefe do Laboratório de Sísmica do Instituto de Geologia e Geofísica, Academia Russa de Ciências, Akademgorodok, Novosibirsk. Atualmente é professor-visitante do LENEP/UENF. Áreas de interesse: processamento de dados sísmicos, problemas inversos e diretos de geofísica, caracterização de reservatórios.

Viatcheslav Ivanovich Priimenko é graduado em Matemática Pura e Aplicada pela Universidade Federal de Novosibirsk (UFN), Rússia, em 1976. Obteve seu doutorado em Física-Matemática em 1990 na UFN. Obteve seu título de Livre Docente em 1997 na UFN. Atualmente é o chefe do LENEP/UENF. Áreas de interesse: problemas diretos e inversos de geofísica e engenharia de petróleo, modelagem numérica, migração e ensino nas áreas de matemática, geofísica e engenharia de petróleo.

Roseane Marchezi Misságia é engenheira civil pela Universidade Católica de Minas Gerais - PUC, Belo Horizonte, Brasil, em 1985. Em 1988 e 2003, obteve os títulos de mestre e doutora em Engenharia de Reservatório e de Exploração, na área de geofísica aplicada, pelo LENEP/UENF. Atualmente, é professora associada do setor de geofísica, no LENEP/UENF. Áreas de interesse: processamento de dados sísmicos, caracterização de propriedades físicas e mecânicas de rochas.

Luis Henrique Amaral é geólogo pela Universidade de São Paulo, São Paulo, Brasil, em 1982. Trabalha na Petrobras desde 1983, na área de aquisição e processamento sísmico. Obteve o título de mestre em 2001 pela Colorado School of Mines, Colorado, EUA, em Geofísica Aplicada à Caracterização de Reservatórios. Atualmente é gerente de processamento geofísico da Petrobras no Rio de Janeiro. Suas áreas de interesse são aquisição e processamento sísmico multicomponente 4D.

• ALEKSEEV AS. 1967. Inverse dynamic seismic problems. In: Some Methods and Algorithms for Interpretation of Geophysical Data. Nauka, Moscow, 9-84.
• ALEKSEEV AS, BESSONOVA EN & MATVEEVA NI. 1979. Inverse Kinematic Problems of Explosive Seismology. Nauka, Moscow, 232 pp.
• BLAGOVESHCHENSKII AS. 2001. Inverse Problems of Wave Processes. VSP Publishing, Utrecht, The Netherlands, 137 pp.
• CERVENÝ V. 2001. Seismic Ray Theory. Cambridge University Press, Cambridge, 713 pp.
• GILL PE, MURRAY W & WRIGHT MH. 1981. Practical Optimization.Academic Press, New York, 465 pp.
• HASKEL NA. 1953. The dispersion of surface waves in multilayered media. Bull. Seismol. Soc. Am., 43: 17-34.
• HIMMELBLAU D. 1972. Applied Nonlinear Programming. McGraw-Hill Book Company, New York, 498 pp.
• JURADO F, CUERT M & RICHARD V. 1995. 1-D layered media: Part 2, Layer-based waveform inversion. Geophysics, 60: 1857-1869.
• KABANIKHIN SI & LORENZI A. 1999. Identification Problems of Wave Phenomena: Theory and Numerics. VSP, Utrecht, The Netherlands, 342 pp.
• KARCHEVSKY AL. 2005. A method for numerical solution of the elasticity system for horizontal layer anisotropic media. Russian Geology and Geophysics, 46: 339-351.
• KENDALL M & STUART A. 1973. The Advanced Theory of Statistics. Vol.2: Inference and Relationship. Charles Griffin & Co. Ltd., London, 723 pp.
• LAVRENTIEV MM, AVDEEV AV, LAVRENTIEV MM Jr & PRIIMENKO VI. 2003. Inverse Problems of Mathematical Physics. VSP, Utrecht, The Netherlands, 275 pp.
• MITROFANOV GM. 1980. Effective representation of the wave field in seismic exploration. Geology and Geophysics (Soviet), 21: 46-58.
• MOLOTKOV LA. 1984. The Matrix Method in the Theory of Wave Propagation in Layered Elastic and Fluid Media. Nauka, Leningrad, 359 pp.
• NOCEDAL J & WRIGHT SJ. 2006. Numerical Optimization. Springer,NY, 664 pp.
• NOLET G (Ed.). 1987. Seismic Tomography, with Applications in Global Seismology and Exploration Geophysics. D. Reidel Publishing Company, Dordrecht, Holland, 386 pp.
• ROMANOV VG. 2002. Investigation Methods for Inverse Problems. VSP, Utrecht, The Netherlands, 280 pp.
• THOMSON WT. 1950. Computation of elastic waves through stratified solid medium. J. Appl. Phys., 21: 89-93.
• URSIN B. 1983. Review of elastic and electromagnetic waves propagation in horizontally layered media. Geophysics, 48: 1063-1081.

Datas de Publicação

• Publicação nesta coleção
  25 Maio 2010
• Data do Fascículo
  Dez 2009

Histórico

• Aceito
  22 Jan 2010
• Recebido
  18 Jun 2008