Utilização das transformadas de Laplace e Fourier-Bessel na modelagem de meios elásticos delgados

Utilização das transformadas de Laplace e Fourier-Bessel na modelagem de meios elásticos delgados

Georgy Mitrofanov^I; Viatcheslav Priimenko^II; Roseane Misságia^III; Luis Amaral^IV

^IInstitute of Geology and Geophysics, Siberian Branch of Russian Academy of Sciences, pr. Koptyuga, 4, Akademgorodok, 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, Imboacica, Macaé, RJ, Brasil. Tel.: [phone]; Fax: [phone] - E-mail: [e-mail]

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

^IVPetrobras / 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

Existem diferentes aspectos que devem ser considerados na modelagem numérica dos campos de onda conectados com reservatórios delgados em subsuperfície. Este trabalho descreve os processos de propagação das ondas sísmicas a partir do sistema das equações de Lamé. Os resultados obtidos são ascomponentes horizontal e vertical do deslocamento, importantes para a análise do levantamento sísmico multicomponente. A fim de proporcionar resultados mais próximos do experimento real, na formulação matemática do problema usamos uma fonte enterrada, do tipo centro de expansão. A solução do problema, construído no domínio espectral a partir das transformadas de Laplace temporal e de Fourier-Bessel espacial, é estudada em detalhes para provar sua aplicabilidade na solução de problemas dinâmicos diretos e inversos da sísmica. Neste caso uma atenção especial é dada à análise da influência da parte real do parâmetro de Laplace nasolução do problema direto, utilizada para a construção de sismogramas multicomponentes sintéticos. Neste artigo, também, discutimos os aspectos computacionais do esquema proposto e realizamos várias análises dos campos das ondas registrados. Os resultados alcançados demonstraram o potencial do método para auxiliarno processo de caracterização de reservatórios delgados em comparação com os métodos potenciais, de raio e de diferenças finitas.

Palavras-chave: sistema de Lamé, reservatório delgado, sistema de observação multi-onda e multicomponente, transformadas de Laplace e de Fourier-Bessel, domínio espectral, modelagem numérica.

ABSTRACT

There are considered various aspects of the numerical simulation of seismic waves connected with thin layer reservoirs. The Lamé system is used for the description of the seismic waves propagation. As a result we obtain both the horizontal and vertical displacement components, which are important for the analysis of seismic fields in the 3D seismic survey. In order to provide results nearer of the real experiment, in the mathematical formulation of the problem we use the deepened source of the type of the center of expansion. The solution of the problem, constructed in the spectral domain using the temporal Laplace and spatial Fourier-Bessel transforms, is studied in details to prove its applicability to the solution of direct and inverse dynamic seismic problems. Special attention is done to analysis of the influence of the real part of the Laplace parameter on multi component seismograms corresponded to the horizontal and vertical components of the elastic field. Also, we discuss the computational aspects of the proposed scheme and carry out several analyses of the registered wave fields. The results obtained demonstrated the potential of the method to help in the process of characterization of thin reservoirs in comparison with other methods: method of potentials, ray tracing method,and finite difference method.

Keywords: Lamé system, thin layer reservoir, multi wave and multi component observation, temporal Laplace and spatial Fourier-Bessel transforms, spectral domain, numerical modeling.

INTRODUÇÃO

Modelagem numérica de campos da onda é uma parte importante do processamento e interpretação de dados sísmicos.A maioria dos métodos de modelagem faz uso da equação da onda acústica. Isto pode estar relacionado a dois fatos. Primeiro, a solução matemática deste problema é bem conhecida e não causa grandes complexidades. Segundo, até recentemente, o interesse da exploração sísmica era isolar e interpretar as ondas P refletidas.

Contudo, hoje, existe um crescente interesse na análise das diversas componentes do deslocamento e métodos de processamento de dados, direcionados para identificação e interpretação das ondas transversais, em particular, ondas convertidas PS. Mas, isto exige o desenvolvimento de métodos matemáticos mais complexos. Neste caso, o uso das equações acústicas para análise da propagação das ondas convertidas é ineficiente, e porisso, para a análise da influência das ondas convertidas é preciso utilizar as equações de Lamé completas.

Neste trabalho, considera-se a solução do sistema de Lamé para modelos horizontalmente estratificados. A fonte sísmica considerada é do tipo pontual de centro de expansão, pois é a que mais se aproxima das fontes explosivas ou canhão de ar comprimido. É importante ressaltar que para este tipo de fonte o sistema de Lamé não pode ser decomposto em equações separadas, e por isto exige uma solução em comum.

Com este obtido, vários trabalhos têm sido direcionados à solução do sistema de Lamé. Os métodos mais difundidos para solução dos problemas da exploração sísmica são os métodos assintóticos (Petrashen, 1978; Petrashen et al., 1982; Cervený, 2001; Chapman, 2004) e de diferenças finitas (Fatianov & Mikhailenko, 1988; Fatianov, 1990; Kelly & Marfurt, 1990). Neste trabalho, usamos os métodos e as idéias formulados em Thomson (1950) e Haskell (1953), veja também Dmitriev (1968); Akkuratov & Dmitriev (1979) e Ursin (1983). A idéia principal é reduzir a solução geral do sistema de equações diferenciais da teoria da elasticidade à solução do sistema de equações do tipo de Riccati. Este tipo de equações apresenta solução na forma analítica para coeficientes constantes por partes, e não acumula erro numérico. A redução é realizada utilizando as transformadas de Fourier-Bessel, com respeito à variável espacial, e Laplace, com respeito à variável temporal. Neste momento é importante que seja efetuado uma análise criteriosa da influência da parte real do parâmetro deLaplace na solução do problema, utilizada para a construção de sismogramas multicomponentes sintéticos.

O algoritmo usado considera inteiramente a especificação horizontal da estratificação do modelo. Por exemplo, a componente relacionada com a superfície de registro pode ser separada de forma muito simples. Isto permite suprimir as ondas diretas e de superfície da solução do problema. Porém, em sua estrutura são preservadas todas as peculiaridades do tipo de fonte adotado.

Uma alternativa de construção da solução completa no domínio espectral é a abordagem baseada no método dos potenciais (Aki & Richards, 2002). Por esta abordagem é possível obter as características espectrais dos diversos tipos de sinaisrefletidos, correspondentes aos modelos elásticos. Apesar desta aproximação não fornecer a solução completa do problema, baseado nela é possível construir com elevada precisão as características espectrais do bloco de camadas delgadas. A vantagem deste procedimento, sobre a solução completa, consiste na proximidade destas soluções com os dados de campo, quando no processamento os diversos sinais refletidos são precisamenteseparados e examinados.

Neste trabalho, fazemos uso do procedimento matemático para solução de problemas sísmicos dinâmicos desenvolvido em Petrashen (1978) e Petrashen et al. (1982), que propõe um algoritmo para solução do problema direto aplicado à solução do problema inverso, considerando uma camada delgada. Após este, um segundo algoritmo mostra como esta solução pode ser estendida para um pacote com poucas camadas delgadas.

A solução do problema direto dinâmico de propagação das ondas sísmicas aplicada ao problema inverso é uma ferramenta importante, e orientada à estimativa dos parâmetros do meio. Neste ponto, investigamos as características da solução do problema direto no domínio espectral, e, também, analisamos em quais casos e o quanto a solução linearizada se aproxima da solução exata. Junto a isto, estudamos a utilização destas características na solução do problema inverso. Além disso, comparamos os resultados dos sismogramas calculados através dos algoritmos propostos e os construídos através dos métodos do raio e de diferenças finitas. Os resultados desta análise comparativa permitiram comprovar a precisão da solução construída e avaliar sua importância no contexto das metodologiasmais divulgadas na sísmica.

FORMULAÇÃO DO PROBLEMA

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é λ, µ, e densidade ρ - 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:

ρ_*π grad_{x, y, z}δ(x, y, z - z_*)g(t),

onde z_* (z_* > 0, z_*≠ 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 ρ_* - 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:

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

A transição 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 ν é a freqüência espacial da transformada de Fourier-Bessel, e p é o parâmetro da transformada de Laplace; p = -α + iw e w = 2πf, f é 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 representados no domínio espectral como:

No problema direto, considera-se a definição das funções (z, ν,p), (z, ν,p), que no domínio z > 0 satisfazem as Eqs. (3) e as seguintes condições de contorno:

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

onde o símbolo [f]_z = f(z + 0) - f(z - 0) representa o salto da função f no ponto z.

As condições definidas pelas Eqs. (4) e (5) são típicas da formulação de problemas matemáticos aplicados à exploração sísmica e garantem a ausência de tensões na superfície livre, decaimento dos deslocamentos quando z → ∞, e continuidade das tensões e deslocamentos nas interfaces das camadas.

Construção da solução completa no domínio espectral

As Eqs. (3) podem ser representadas na forma matricial

onde

O símbolo * significa a transposição da matriz, e M, N, K são definidas pelas seguintes fórmulas:

onde

e c é o tipo da onda utilizada: p - para a onda primária e s - para a onda secundária;

Representações similares são usadas para as condições de contorno

e para as condições nas superfícies de descontinuidade z = z_k, k = 1, 2,..., n:

Agora, consideremos o sistema auxiliar:

e o vetor da função a priori G,

G ≠ 0.

A solução do problema U(z, ν, p), apresentada nas Eqs. (7)-(9), é construída utilizando o método de potenciais, que consiste na representação do campo vetorial através da soma de dois campos vetoriais - o potencial e o rotacional. O campo potencial possui o potencial escalar (z, ν, p), e o campo rotacional - o potencial vetorial (z, ν, p). Os potenciais são definidos através das seguintes equações diferenciais:

As funções , estão relacionadas aos espectros das funções u, w através das seguintes fórmulas:

Introduzimos a matriz S = (s_ij ,1 < i, j < 2), que estabelece uma relação entre os potenciais e suas primeiras derivadas:

Substituindo as Eqs. (12) em (10) é possível constatar que os elementos desta matriz obedecem ao sistema de equações diferenciais ordinárias não-lineares de primeira ordem do tipo de Riccati:

Quando r_p, r_s são constantes ou constantes por partes, a solução pode ser representada na forma explícita ou por fórmulas recursivas, respectivamente. Ao mesmo tempo, através das duas últimas equações é possível concluir que s₁₂: s₂₁ = const ou, mais exatamente, s₁₂ = ν²s₂₁ . As funções v_p(z), v_s(z) e ρ(z) são constantes por partes, com respeito à variável z ∈ (0, ∞); r_p(z), r_s(z) apresentam as mesmas propriedades. Em virtude disto, as Eqs. (13), em cada intervalo

z ∈ (z_k–1, z_k), k = 1, 2,..., n,

têm uma solução representada na forma explícita. Devido às condições de contorno definidas pelas Eqs. (7), para z ∈ (z_n, ∞), a matriz S têm forma diagonal:

Portanto os elementos da matriz s_ij, i, j = 1, 2, são conhecidos em z = z_n + 0. As Eqs. (11) e (13) permitem obter as condições de transmissão para as funções s_ij e definir s_ij(z_n - 0), respectivamente. A seguir, utilizando os últimos valores como condição de fronteira, é possível determinar s_ij(z_n–1 + 0). Apartir das condições de transmissão calcula-se s_ij(z_n–1 - 0). Repetindo este procedimento camada por camada, calculamos os valores s_ij(0). Agora, utilizando as Eqs. (7), (11) e (12), é possível obter:

onde

A solução do problema auxiliar expresso nas Eqs. (7)-(9) pode ser generalizada da seguinte forma:

1. Cálculo de s_ij, i, j = 1, 2, na superfície z = 0, de baixo para cima.

2. Cálculo de valores dos potenciais na superfície z = 0.

3. Cálculo de valores dos espectros.

Agora, consideremos a solução do problema principal apresentado nas Eqs. (3)-(5). Nesse caso

caracteriza uma fonte externa. Para z₀ Î (z_k-1, z_k), k = 1, 2,..., n - 1, vamos buscar a solução do problema da seguinte forma:

onde

Cada função U_k, k = 1, 3, 4, deve satisfazer as seguintes equações:

Substituindo a Eq. (15) na Eq. (6), utilizando as Eqs. (17) e (18), igualando os coeficientes com

obtemos as seguintes igualdades:

Para cada função vetorial U_k, k = 1, 3, 4, é preciso definir os potenciais escalar e vetorial _k, _k, k = 1, 3, 4, para os quais, utilizando relações similares às Eqs. (12), possibilita introduzir as matrizes S, X, Y, correspondentes. Usando igualdades similares às Eqs. (11), e a partir da Eq. (15), obtemos:

onde x_ij, y_ij, i, j = 1, 2, são os elementos das matrizes X, Y. Introduzimos a função vetorial

Através de relações similares às Eqs. (12), temos:

e da Eq. (7), onde G = 0, segue-se que:

B_S Φ₁(0, ν, p) = B_Y Φ₄(0, ν, p).

Então das condições similares às Eqs. (11) e (12) e da Eq. (15) obtemos

onde

Aqui j_km, k, m = 1, 2, representam os elementos damatriz J = S, Y.

As fórmulas obtidas estão corretas até a determinação da matriz Y, em z = z₀ . Definimos

(se o meio abaixo de z = z₀ é um semi-espaço compreendendo os valores dos parâmetros da camada m, então a matriz Y corresponde a uma solução particular U do problema original, e U → 0 quando z → ∞). Assim, o valor U em z = 0 é completamente definido. Observe que S ≡ X quando z ∈ (z₀, ∞). Agora, supondo que z₀∈ (0, z₁) é possível mostrar que Y(z, ν, p) = Y(z₀, ν, p), z ∈ (0, z₁), e a matriz P pode ser representada da seguinte forma:

Feito isso, está concluída a construção da solução correspondente ao problema direto. Como resultado, quantificamos as componentes do vetor U, na forma de espectros bidimensionais (z, ν, p), (z, ν, p).

A partir deste método é possível propor um algoritmo para estimar, com razoável precisão, os valores das funções, considerando um grande número de camadas (não menos de 50). O algoritmo proposto pode ser usado para a simulação de dados sísmicos de acordo com o seguinte diagrama:

1. aplicação de transformadas de Fourier-Bessel e de Laplace em variáveis apropriadas as equações e as condições iniciais e de contorno;

2. solução das equações correspondentes no domínio espectral;

3. realização das transformadas inversas em função das variáveis ν e p.

O passo seguinte é fornecer os dados de entrada, tais como: descrição do modelo, dos parâmetros do arranjo, profundidade da fonte, e da forma do sinal. No final desta simulação obtemos dois sismogramas correspondentes às componentes horizontal e vertical u(z, r, t), w(z, r, t).

Características espectrais de reflexão

A solução obtida permite a construção explícita das características espectrais de vários modos de ondas refletidas em alvos com estruturas horizontalmente estratificadas. Do ponto de vista prático, a expectativa é obter as características das ondas PP e PS, construídas para o caso de incidência de uma onda planalongitudinal num pacote de camadas delgadas.

As construções apresentadas demonstram o fundamento para estimar os parâmetros da região de camadas delgadas com base no método clássico de análise AVO. A metodologia proposta se diferencia do método clássico de análise AVO, visto que na solução do problema inverso não são consideradas as reflexões do topo e da base da camada delgada, impossíveis de serem separadas pelo campo de onda, e sim, a reflexão da camada delgada como um só objeto. É importante ressaltar que a substituição da interface entre os dois semi-espaços elásticos por uma camada delgada, e depois, por um conjunto de camadas delgadas, complica consideravelmente a solução do problema inverso. Neste caso, não estamos interessados na divisão das reflexões no campo de onda. Nesta metodologia, usamos as características espectrais do pacote de camadas delgadas, em vez dos coeficientes de amplitude da reflexão. Por isso, é necessário efetuarmos o estudo dos problemas ligados à interface entre dois semi-espaços (ou camadas espessas), não como um contato rígido, mas como um objeto mais complexo, representado por camada delgada (ou conjunto de camadas delgadas). A importância disto para aplicações práticas está relacionada ao fato de que a maioria dos reservatórios de óleo e gás serem formados por pacotes de camadas delgadas, cujo estudo, continua sendo um problema básico da prospecção sísmica.

O potencial da onda longitudinal plana incidente pode ser representado da seguinte forma

onde

k é o número de onda e é o ângulo de incidência.

Para a construção da solução correspondente usamos a representação da onda plana na forma de uma série de Fourier com respeito às funções de Bessel, veja Mittra & Lee (1971). A partir desta representação e do método proposto, podemos construir os potenciais das ondas PP e PS refletidas e refratadas em uma camada delgada, e obter os coeficientes complexos de reflexão, que define as características espectrais de cada tipo de onda. A seguir, representamos os coeficientes de reflexão/ refração, como uma superposição de ondas planas, refletidas no topo e na base da camada. Neste caso, a amplitude das ondas que sofreram mais que uma conversão torna-se desprezível. São fracas, também, as reverberações das ondas provenientes de mais de uma reflexão na interface com baixo contraste de impedância elástico. Por conta disto, vamos considerar somente os coeficientes de reflexão/refração da primeira interação das ondas com o meio. Então, usando as fórmulas de Zoeppritz (Zoeppritz, 1919) linearizadas para uma interface, definimos as características espectrais e da camada delgada, correspondentes aos tipos de ondas indicados, veja Figura 1.

onde

e

Para um modelo composto de diversas camadas, a solução pode ser obtida de forma análoga:

onde N é o número de camadas no pacote, e as expressões para e podem ser escritas na seguinte forma:

Tais fórmulas serão usadas nos procedimentos da simulação e solução do problema inverso de objetos delgados. Para os expoentes que tomam parte destas fórmulas usamos as seguintes representações:

Entretanto, a simplicidade do esquema proposto não evitou os problemas. Sendo assim, para garantir uma realização prática eficaz é necessário examinar as características especiais da solução construída no domínio espectral, e as características especiais das transformações adotadas. A seguir, com base no algoritmo proposto, analisamos as características especiais e as correspondentes transformadas, e apresentamos os sismogramas calculados para os modelos de uma camada delgada e, também, do bloco de camadas delgadas.

ANÁLISE DA SOLUÇÃO NO DOMÍNIO ESPECTRAL

Com base no algoritmo citado anteriormente, foram criados vários procedimentos. Estes permitiram calcular efetivamente, em termos de rapidez e precisão, as características espectrais bidimensionais (z, v, p), (z, v, p), correspondentes às componentes horizontal e vertical do vetor de deslocamento. Para implementação deste algoritmo, é necessário que seja fornecido os seguintes dados: descrição do modelo, parâmetros de espaçamento, profundidade da fonte e forma do sinal gerado.O produto resultante é a criação de um programa, que permite obter, com base nos espectros, os dois sismogramas correspondentes às componentes horizontal e vertical: u(z, r, t), w(z, r, t).

Um número razoável de modelos foi usado para testar e validar as propriedades dos algoritmos construídos. A Figura 2(a) representa um modelo, cuja camada delgada corresponde a um arenito saturado por gás, posicionada a 500 metros de profundidade. Apesar da simplicidade, este modelo permite investigar detalhadamente os procedimentos de cálculo da solução do problema inverso sob vários aspectos. Em primeiro lugar, no domínio do tempo, permite decompor com exatidão o caráter da reflexão e estimar a precisão do cálculo de cada componente do campo de onda. Em segundo lugar, no domínio espectral, lidamos com a análise da interferência para identificar a solução construída. Como demonstra este trabalho, a complicação do modelo não eleva a precisão das representações. Além disso, este modelo permite formular exigências em relaçãoà transformada inversa, que indicam o nível de precisão dos sismogramas construídos.

Na realização das simulações utilizamos o impulso de Ricker,com freqüência dominante de 30 Hz e 50 Hz. A freqüência de 30 Hz é uma das mais usadas nas geometrias de aquisição dedados reais. A freqüência de 50 Hz foi escolhida para testar a capacidade de resolução do campo de onda dos algoritmos. Para cálculo dos sismogramas e verificação dos algoritmospara solução do problema inverso, usamos 48 receptores igualmente espaçados de 25 m. A posição do primeiro receptor coincide com a posição da fonte ou pode estar deslocada de 25 m em relação à posição da fonte. Sendo assim, o comprimento do offset alcança 1200 m. A fonte está localizada a 15 m de profundidade da superfície livre ou na mesma.

A seguir, descrevemos as principais modelagens realizadas. Inicialmente, são ilustradas as particularidades dos espectros obtidos para o modelo mais simples. Na Figura 3 é apresentada a parte real do espectro da componente vertical (p, ν) ≡ (0, p, ν), usando α = 0.01. Estas figuras mostram detalhes do comportamento dos espectros 2D, obtidos na solução do problema direto. O espectro da componente horizontal (p, ν) ≡ (0, p, ν) contém uma estrutura similar às encontradas no espectro da componente vertical. Através da análise da Figura 4, é possível visualizar que em todos os espectros, calculados para os outros modelos, observamos a mesma situação. Deste modo, a fim de investigar maiores detalhes no cálculo dos espectros, adotamos o intervalo regular de 0 a 0.5 para a freqüência espacial ν, e diferentes intervalos de freqüência temporal ω.

Assim, na Figura 3(a) estão representados os espectros para o intervalo da freqüência temporal de 1 até 151 Hz, com passo de 3 Hz. Nesta figura é possível observar o caráter irregular dos espectros, correspondente ao intervalo de freqüências analisados. Diante disto, a estrutura geral do espectro construído pode ser separada em zonas de gradual a bruscas variações de valores do espectro. Uma análise destas zonas demonstra que estas bruscas variações são bem localizadas e deslocam-se em função da freqüência espacial, respondendo ao incremento da freqüência temporal como se fosse criada uma zona de expansão das singularidades do espectro 2D. Sendo assim, para valores fixos da freqüência espacial ν ou temporal ω, são definidos intervalos, onde os espectros têm caráter gradual. Intervalos mais restritos de freqüência temporal, Figura 3(b), permitem observar a diferença no caráter das variações dos valores do espectro em função da freqüência espacial ν. A partir disto, separamos três zonas, caracterizadas pelas diferentes variações do espectro. Na primeira zona, as variações da freqüência temporal referem-se aos valores de freqüência espacial de 0.1 a 0.14, e correspondem ao caráter descontínuo do espectro, acima citado. A segunda zona corresponde aos valores ν de 0.24 a 0.29 e apresenta caráter claro de ruídos, indicando certa instabilidadenas soluções construídas.

Esta instabilidade tem caráter randômico e variações desprezíveis, se comparada com a variação total do espectro. Neste caso, o parâmetro α pode atuar como um parâmetro regulador. Finalmente, a terceira zona corresponde aos valores ν de 0.3 a 0.4. As singularidades do espetro observado apresentam mais explicitamente a descontinuidade da estrutura. Através desta análise investigamos a natureza das singularidades, e mostramos sua conexão com a parte da solução determinada pela presença da superfície livre (ondas primárias e de superfície). A exclusão das componentes relacionadas à superfície de observação possibilitou remover completamente as singularidades correspondentes aos espetros bidimensionais. Isto é feito subtraindo da solução completa a parte conectada com a superfície livre. A Figura 3(c) demonstra o comportamento dos espectros nas principais zonas com fortes variações. Na Figura 3(b), a escala ampliada do eixo mostra os valores do espectro. É visível, quenestas zonas, a estrutura do espectro se aproxima às descontinuidades ou δ-funções. Uma hipótese é supor que estas singularidades estão conectadas com as diferentes soluções do sistema de Lamé.

A Figura 4 comprova a similaridade entre a estrutura dosespectros 2D obtidos em diferentes modelos. Das imagens dos três espectros citados, correspondentes aos modelos analisados, é possível notar que a estrutura na região com propriedades de descontinuidade do espectro é pouco alterada. Mas, o comportamento bem localizado desta estrutura, permite combinar as freqüências espacial e temporal, onde estas singularidades estão ausentes. Então, a partir da combinação de freqüências é possível construir espectros bem suavizados. Tal procedimento deve ser usado na solução do problema dinâmico inverso por garantir a alta convergência.

Os resultados da Figura 5 ilustram mais uma característica importante da solução construída no domínio do espectro - a dependência da estrutura dos espectros do parâmetro α. A comparação das Figuras 4(a) e 5(a) demonstra, mediante a variação do parâmetro α de 0.01 a 1, o quanto pode ser essencial a alteração do espectro para a solução geral. A Figura 5(b) ilustra que as 10 alterações realizadas na solução truncada, de fato excluem todas as singularidades, e suavizam o espectro das funções. Semelhantes simulações, efetuadas em diferentes modelos, permitiram tecer considerações. Primeiro, o parâmetro α define uma função amortecida com respeito de variável t, e a forte influência nos sinais sísmicos, associados com os maiores tempos de trânsito, indica a possível naturezadas descontinuidades do espectro das soluções adequadas às respectivas ondas, por exemplo, reverberações, etc. Segundo, utilizando o algoritmo proposto, e através da aplicação do parâmetro α, fica comprovado que a suavização dos espectros pode ser realizada efetivamente no processo de construção de funcionais objetivos, utilizados na solução de problemas inversos dinâmicos da sísmica (Priimenko et al., 2007). Terceiro, a suavização do espectro e eliminação das descontinuidades indicam a possibilidade de aplicação do parâmetro α como um parâmetro de regularização no acoplamento das diferentes soluções, e, também, dos espectros, obtidos por meio de dados reais. Além disto, observamos a possibilidade de estimar oserros presentes na solução do problema, através da substituição das fórmulas exatas pelas aproximada. Exemplos destas estimativas, utilizando as fórmulas aproximadas e os potenciais simples, são mostrados na Figura 6. Note nas Figuras 6(a) e 6(c) a presença das características espectrais construídas para asondas longitudinais e convertidas do modelo da Figura 2(a).

Também, apresentamos os resultados da análise comparativa do percentual de erro entre a solução exata e a aproximada. A efetividade desta análise é assegurada quando não representamos a solução exata do campo de onda total, mas somente para algumas de suas componentes. Na análise de dados reais, este procedimento é de suma importância para separar com exatidão estes componentes. Um exemplo disto são as ondas longitudinais e convertidas monotípicas. A investigação realizada para estimar os erros em função da freqüência (ou a proporção H/λ) e do ângulo de incidência das ondas incidentes no refletor (ou o ponto médio do afastamento entre fonte e receptor X/2) em diferentes modelos, permite que sejam tecidas algumas conclusões. Primeiro, as alterações das características espectrais, provenientes das fórmulas aproximadas, não superam 2%. Segundo, a comparação entre os erros de aproximação dos valores exatos e os aproximados tem caráter periódico em função da freqüência, e tendem a aumentar com o afastamento da fonte. Terceiro, existem zonas de freqüências, onde os erros de tais aproximações são mínimos. Junto a isto, nas ondas longitudinais a minimização dos erros, considerando o caso de substituição das características espectrais exatas pelas aproximadas, ocorre nas condições H/λ ≈ 1/4, 3/4 e etc.; nas ondas convertidas PS esta condição não é válida. Isto ocorre como conseqüência da grande influência que os coeficientes de reflexão exercem sobre as ondas longitudinais, assim como nas transversais.

ANÁLISE DA SOLUÇÃO NO DOMÍNIO TEMPORAL

Do ponto de vista sísmico, a solução do problema inverso dinâmico no domínio espectral dificilmente é submetida à análise. Por esta razão, grandes esforços foram feitos para o desenvolvimento de procedimentos estáveis, que assegurem a qualidade dos sismogramas nos variados domínios (espectral - tempo e espaço). Como resultado, é possível mostrar que os sismogramas, construídos através da aplicação do procedimento proposto, oferecem uma solução bem próxima da ideal, mesmo para o caso da modelagem de reservatórios com espessuras subsísmicas. Tal procedimento elevou o grau de confiabilidade da solução do problema inverso no domínio espectral.

Partimos do algoritmo usado na transição do domínio temporal para o espectral. Observamos os princípios gerais do algoritmo e assinalamos as mais importantes características. Conforme a aproximação inicial do problema, para obtenção das componentes u(z, r, t), w(z, r, t) é preciso realizar duas transformações inversas, uma em função da freqüência espacial e outra temporal, utilizando os valores dos espectros calculados. Uma breve análise das respectivas transformadas indica que ao realizarmos a transformada inversa, em função do espaço, é preciso calcular duas integrais com as funções de Bessel, o nulo e o primeiro, que respondem às correspondentes transformadas diretas para as componentes horizontal e vertical. Ao realizar a transformada inversa em função da freqüência temporal, que responde à transformada de Laplace, é possível utilizar a transformada geral de Fourier. O desmembramento das transformadas permitiu criar um algoritmo mais efetivo. Assim, no cálculo das integrais com as funções de Bessel foram aplicadas fórmulas de quadratura de Gauss com 64 ordenadas, assegurando uma boa aproximação para os espectros fortemente irregulares, veja, por exemplo, (Demidovich & Maron, 1987). Na construção da transformada inversa de Fourier, utilizamos o método de Failon para calcular a respectiva integral trigonométrica. Este tipo de integral apresenta alto nível de precisão, veja (Krylov, 2006), e apesar de levar em consideração a homogeneidade dos dados de entrada (traços ou espectros), não é menos rápida que a FFT (transformada rápida de Fourier). A percepção da estrutura dos espectros 2D é um ponto de extrema relevância para aplicação dos procedimentos dos cálculos indicados. Tal análise possibilita construir de forma correta o algoritmo de cálculo das referidas integrais. Além disso, também, permitiu escolher corretamente os parâmetros de discretização dos espectros e aplicar efetivos filtros de suavização, estabilizando as transformadas inversas. Sendo assim, após uma análise dos espectros 2D é possível escolher o algoritmo de cálculo das integrais, correspondentes à transformada de Fourier-Bessel, com respeito à freqüência espacial, e determinar o respectivo passo de discretização. Por exemplo, partindo dos procedimentos formais, semelhante ao caso da transformada de Fourier, podemos estimar o requerido passo com respeito à freqüência espacial ν. Se precisarmos calcular os traços sísmicos para posicionamento dos receptores com abertura máxima de L = 2500 m, o suposto passo Δν deve estar na ordem de 1/L = 0.004. Ao mesmo tempo, conforme a Figura 5(c), onde se tinha passo semelhante, nos domínios com grande variação dos espectros a grandeza do passo Δν deve ser bem reduzida para garantir o cálculo correto das integrais. As simulações realizadas mostraram que nos domínios indicados é preciso diminuir em quase 100 vezes o valor do passo. Considerando as fórmulas de quadratura de Gauss, este procedimento assegura a estabilidade durante o cálculo das respectivas integrais. Abaixo, durante a tentativa de construção dos requeridos traços sísmicos serão mostradas as conseqüências do dimensionamento errado do passo Δν. Com respeito à transformada inversa de Fourier os procedimentos formais não garantem um bom resultado. Por exemplo, analisando a discretizaçãonos domínios tempo-freqüência para aplicação da transformada de Fourier, é fácil produzir uma estimativa aproximada para os respectivos passos e intervalos de discretização. Assim, se quisermos obter o registro do modelo até 2.5 s com intervalo de amostragem de 2 ms, é necessário adotar o intervalo de freqüência até 250 Hz, cujo intervalo de amostragem é de 0.4 Hz. Ao mesmo tempo, uma análise realizada nos espectros 2D, assim como em 1D, obtidos ao final da integração com respeito a ν (veja Fig. 8), mostra um espectro mais irregular em função da freqüência temporal ω do que esperava-se. Como indicado acima, isto está relacionado à presença da informação sobre todos os possíveis sinais na solução construída do problema.Por isso, temos que usar o parâmetro α, para diminuir sua influência, ou incrementar substancialmente a quantidade dos pontos de discretização do intervalo de freqüência.

A consideração de singularidades dos espectros indicadosna escolha de algoritmos efetivos para cálculo das requeridas integrais permite calcular sismogramas para as componenteshorizontal e vertical do campo de onda, correspondentes a modelos de camadas delgadas de qualquer complexidade. A solução do problema direto para modelos com grande número de camadas eleva, essencialmente, o tempo necessário para obtenção da solução. Além disto, este tipo de modelo aumenta em muito a complexidade dos espectros, incrementando o número de intervalos de discretização exigidos para a realização das transformadas inversas de Fourier-Bessel e Laplace e, também, à aplicação dos procedimentos específicos. Ao mesmo tempo, ao implementarmos o algoritmo desenvolvido, é possível obter sismogramas de alta qualidade, que de fato aproximam os dados sintéticos à solução analítica do problema. Quando se trabalha com modelos compostos é bom considerar que a complexidade da interpretação dos sismogramas sintéticos obtidos se aproxima à complexidade da interpretação dos sismogramas reais dos campos de onda. Neste caso a reflexão obtida tem um caráter interferencial muito complexo. Por isso, do ponto de vista de uma verificação detalhada sobre a capacidade do algoritmo, os sismogramas mais importantes são representados pelos resultados que foram obtidos em modelos relativamente simples, com pequeno número de camadas. Este sismograma multicomponente, embora elementar, mostra a possibilidade de identificar os diferentes tipos de componentes das ondas. Por exemplo, através da componente vertical (com certo afastamentodos receptores da fonte), foram observadas as ondas convertida e transversal, além da onda refletida geral. E as componentes horizontal e vertical destes mesmos sismogramas permitiram analisar as particularidades na formação da dinâmica das ondas. Assim, revelam-se as variações dinâmicas das diferentes componentes do campo de onda nos pontos de registro, correspondentes aos diferentes afastamentos da fonte. Por exemplo, analisando a componente horizontal observamos a variação da amplitude e da dinâmica das ondas convertida etransversal, mais afastadas da fonte. Junto a isto, quanto mais distante estiver à onda convertida da fonte, mais lento é o seu amortecimento, em comparação com a onda transversal. Além disso, os sismogramas obtidos permitiram estimar a grandezadas amplitudes das diferentes ondas nas diferentes componentes. O que é considerado como um resultado interessante na aquisição de dados multicomponentes.

As Figuras 7 e 8 mostram os resultados das transformadas inversas para espectros que foram construídos como solução do problema direto em modelo contendo uma só camada (veja a Fig. 2(a)). Com isso, as transformadas inversas, correspondem aos parâmetros usados na Figura 7. Nesta figura, para maior detalhamento, foram mostrados somente cinco traços, situados perto da fonte (o grupo de traços de cima corresponde à componente horizontal, e o de baixo à componente vertical). Observe que as componentes correspondentes à superfície livre foramexcluídas desta solução. Isto permitiu analisar minuciosamente algumas características das ondas refletidas, que podem representar o maior interesse da interpretação, e, também, estimar o nível das componentes randômicas, associadas com o procedimento de cálculo das transformadas inversas aplicadas. A Figura 7 representa um sismograma construído com base na solução do problema no domínio espectral. Esta figura mostra que através da solução do problema no domínio espectral épossível construir traços em regiões onde o nível da instabilidade computacional é desprezível. Esta situação pode ser visualizada nos receptores mais próximos da fonte, diminuindo significativamente em direção aos mais distantes da fonte. Comisso, é possível minimizar a presença de componentes randômicas nas componentes verticais, relacionadas com a estabilidade do algoritmo.

Conforme demonstra os resultados da Figura 7, nos traços obtidos pode ser observados claramente, não somente os diversos modos de ondas, incluindo a onda convertida refletida, mas também as reverberações. Assim, através da componente vertical ω, próximo a 1.1 s, observamos a primeira reverberação da onda longitudinal. Na componente vertical podemos identificar outras ondas (convertida e transversal) que começam a se formar em torno da região com 0.9 s, mesmo diante de pequenos afastamentos. Na componente horizontal nitidamente revelam-se todos os tipos de ondas. Com isso, a amplitude aumenta à medida que se afasta da fonte, o que corresponde à natureza das ondas e ao tipo da fonte explorada (centro de dilatação). Convêm, também, notar, que para a componente horizontal, a relação entre as amplitudes dos diversos tipos de ondas difere totalmenteda componente vertical.

Os sismogramas da Figura 8 apresentam o campo de onda geral que inclui as ondas direta e de superfície. Através da análise das componentes horizontal e vertical é possível identificar os detalhes da formação dos diferentes tipos de ondas. Além disso, pelos sismogramas é possível investigar a variação dinâmica das componentes do campo de onda nos pontos de registro, quecorrespondem aos diferentes afastamentos entre fontes e receptores. Particularmente, através da componente horizontal é possível observar que a amplitude e a dinâmica das ondas convertida e transversal variam com o afastamento da fonte. Junto a isto, a onda convertida é amortecida mais lentamente do que a transversal, fato teórico bem conhecido. Na Figura 8, através da componente vertical podemos observar a onda refletida, a convertida(de 0.9 s a 1.2 s) e a transversal (de 1.3 s até 1.9 s). Ao mesmo tempo, estes sismogramas mostram o quanto a amplitude destas ondas é mais evidente na componente horizontal do que na vertical, em comparação com onda longitudinal refletida. Assim, pela componente horizontal a amplitude das ondas convertidas é comparável com a da onda longitudinal.

O último sismograma serve como mais uma justificativa do uso de instrumentos multicomponentes para a separação das referidas ondas, exigida na solução do problema inverso. Além disso, nos sismogramas apresentados estão bem destacadas, para cada tipo de onda, as regiões com crescimento da amplitude. A estrutura do campo pode ser útil no nível de planejamento que utiliza os receptores multicomponentes. Assim como, durante o processamento dos dados reais, visando à separação dos diferentes tipos de ondas. Para validação dos sismogramas construídos é razoável compará-los com sismogramas calculados a partir de métodos conhecidos. O método de modelagem por diferenças finitas é o mais comum entre os que podem ser considerados. As possibilidades e propriedades deste método sãobem conhecidas. Nas Figuras 8(c) e 8(d) mostramos os resultados dos sismogramas sintéticos obtidos pelo método de diferenças finitas. Note que é possível construir todos os tipos de campos de ondas, mas com qualidade inferior ao do algoritmo proposto. Assim, as características dinâmicas dos sinais, correspondentes às ondas, não são bem reveladas, por força da peculiar "difusão do impulso", sofridas por muitos algoritmos baseados no método de diferenças finitas. Além do método de diferenças finitas, também analisamos o método de traçamento de raio.

No modelo da Figura 2(a) usamos o método de raio para calcular somente a componente vertical. Para calcular ambas componentes, horizontal e vertical, precisaríamos usar fontes com diferentes potências. Os resultados obtidos são apresentados na Figura 9(b). Levando em consideração a impossibilidade de modelar as ondas de superfície dentro dos programas do método de traçamento de raio disponíveis, estes foram excluídos da solução do problema direto no algoritmo proposto. Na Figura 9(a) apresentamos um sismograma correspondente à componente vertical, obtido pelo algoritmo proposto. A comparação entre estes sismogramas mostra que apesar do método de traçamento de raio permitir formar as ondas convertida e transversal com propriedades cinemáticas corretas (tempo de trânsito das ondas), Figura 9(b), as amplitudes não correspondem às obtidas com a solução exata do problema direto, Figura 9(a). Isto pode ser justificado pelo uso de um tipo diferente de fonte. Desta forma, ao aplicar o método de traçamento de raio, temos que utilizar uma fonte que combine as forças horizontal e vertical, e assim obter as ondas convertida e longitudinal. Ao mesmo tempo, no algoritmo proposto para solução do problema direto, usamos uma fonte do tipo do centro de dilatação, posicionada abaixo da superfície livre. No sismograma, construído pelo método de traçamento de raios, é visível que a fonte combinada levou ao surgimento da componente de onda transversal. Este último efeito nunca foi observado em situação real.

A Figura 10 demonstra os resultados mais complicados de aplicação do algoritmo proposto para cálculo dos sismogramas correspondentes às componentes horizontal e vertical.A Figura 2(c) apresenta os parâmetros do modelo multilaminar usado, e da subcamada. Como nos exemplos anteriores, durante os cálculos dos sismogramas foram extraídas as componentesdo campo de onda, correspondentes à extremidade livre (onda direta e a de superfície).

Na execução das transformadas inversas, foram usados os mesmos parâmetros empregados nas modelagens mais simples. Por isso, não é possível garantir uma boa qualidade do campo de onda calculado. Mas aqui não será efetuado este tipo de análise, que se justifica em sismogramas provenientes de dados reais, cuja finalidade é garantir uma melhor interpretação e solução de problemas geológicos concretos. No presente caso, a importância está em estimar parâmetros de objetos complexos através do algoritmo proposto.

Em princípio, os sismogramas obtidos são verossímeis.Neles estão presentes os diferentes tipos das ondas. Todas as ondas possuem uma natureza complexa interferencial, o que está relacionada à complexidade do meio. As principais características do comportamento dos diferentes tipos das ondas nas diversas componentes coincidem com a modelagem mais simples, analisada anteriormente. Assim, através da componentehorizontal podemos observar pequenas amplitudes nas ondas longitudinais. Ao mesmo tempo, na estrutura do campo de onda obtida surgem características adicionais, que o diferem do campo de onda para a camada delgada. Por exemplo, uma redistribuição de energia pode ser observada nas ondas convertida e transversal, comportamento peculiar na componente horizontal. Em uma análise rápida, esta característica pode levar a sério erros, e este deve ser considerado em casos reais. Além disso, para um pacote com considerável número de camadas delgadas, também é possível observar pela componente horizontal reflexões mais significantes do que a obtida em interfaces individuais ou camadas delgadas isoladas. Baseado nisto, pode-se esperar que diante da análise dos reservatórios de camadas delgadas, as componentes horizontais forneçam informações através das ondas convertidas.

A interpretação mais detalhada do caráter da reflexão dentro da presente modelagem, não foi realizada. Em nossa opinião, uma observação detalhada do caráter da reflexão do modelo e sua interpretação devem ser realizadas em conjunto com a solução do problema inverso. É preferível, também, compará-los com situações reais observadas. No último caso, existe a oportunidade de estimar: o quanto o campo modelado corresponde à situação real e quais de suas características têm maior significado.

CONCLUSÕES

O algoritmo proposto neste trabalho tomou como base a transformação do sistema completo de Lamé, no domínio espectral, usando as transformadas de Fourier-Bessel com respeito à variável r, e a de Laplace com respeito à variável t. Isto constitui um esquema eficaz para solução do sistema construído, e permite considerar o algoritmo como uma alternativa ao esquema de diferenças finitas, amplamente utilizadas na modelagem numérica de meios horizontalmente estratificados. Tal procedimento possui vantagens relevantes:

1. Exclusão das componentes que correspondem à superfície de observação, remoção dos efeitos de superposição das ondas de superfície e primária nas ondas refletidas nos objetos alvos.

2. A duração e a qualidade do trabalho não são determinadas pelo intervalo de discretização do modelo do meio; por isso o método proposto demonstra ser eficaz em modelos que representem de forma realística algumas situações estratigráficas encontradas em reservatórios, como: acunhamento, afinamentos estratigráficos de corpos de areia e bordas de estruturas canalizadas, que provocam um efeito significativo na amplitude da onda refletida.

3. Demonstra elevada estabilidade, no que diz respeito à acumulação de erros numéricos e preservação da forma inicial do sinal, possibilitando a obtenção de sismogramas sintéticos mais próximos daqueles obtidos com a solução teórica exata do problema.

Nesta pesquisa demonstramos o potencial do algoritmo proposto para auxiliar no processo de caracterização de reservatórios constituídos por corpos com espessuras subsísmicas. As considerações aqui efetuadas sugerem a continuidade desta pesquisa através da aplicação do algoritmo proposto na análise de imagens das ondas e investigações multicomponentes. Além disso, levando em consideração a elevada precisão dos sismogramas sintéticos apresentados, torna-se interessante propor um algoritmo para solução do problema inverso no domínio espectral, baseado no traçamento de raios de diferentes modos de ondas. Espera-se que o uso deste algoritmo de inversão, nas etapas de exploração e desenvolvimento de um campo, produza imagens mais fidedignas da geologia em subsuperfície. E através da integração de todos os dados disponíveis possa ser gerado modelos mais precisos, otimizando assim, a localização de zonas produtoras e principalmente, aumentando o sucesso na previsão do comportamento de um campo.

AGRADECIMENTOS

Os autores agradecem aos revisores anônimos pelas correções realizadas no manuscrito original. Os autores também agradecem à PETROBRAS pelo apoio financeiro dado para a execução deste trabalho.

Recebido em 18 junho, 2008 / Aceito em 23 junho, 2009

Received on June 18, 2008 / Accepted on June 23, 2009

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 o seu título de Mestre em Matemática na UFN em 1978. 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. Á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 e sísmica multicomponente e 4D.

• AKI K & RICHARDS PG. 2002. Quantitative Seismology. 2^nd edition. University Science Books, Sausalito, California. 596 pp.
• AKKURATOV GV & DMITRIEV VI. 1979. A method for calculation of stationary elastic oscillations in stratified media. In: Numerical Methods in Geophysics. Moscow State University Publ., Moscow, Russia, 3-12.
• BAN-MENAHEM A & SINGH SJ. 1981. Seismic Waves and Sources. Springer-Verlag, New York. 1108 pp.
• ERVENÝ V. 2001. Seismic Ray Theory. Cambridge, New York: Cambridge University Press. 713 pp.
• CHAPMAN CH. 2004. Fundamentals of Seismic Wave Propagation. Cambridge: Cambridge University Press. 646 pp.
• DEMIDOVICH BP & MARON IA. 1987. Computational Mathematics. Mir Publishers, Moscow. 688 pp.
• DMITRIEV VI. 1968. General method for calculus of electromagnetic field in stratified medium. In: Numerical Methods and Programming. Moscow State University Publ., Moscow, Russia, 55-65.
• FATIANOV AG. 1990. Semi-analytical method for the solution of direct dynamical problems in stratified media. Dokl. AN SSSR, 310(2):323-327.
• FATIANOV AG & MIKHAILENKO BG. 1988. A Method for calculation of non-stationary wave fields in non-elastic stratified-inhomogeneous media. Dokl. AN SSSR, 301(4): 834-839.
• HASKELL NA. 1953. The dispersion of surface waves on multilayered media. Bulletin of the Seismological Society of America, 43: 17-34.
• KELLY KR & MARFURT KJ (Eds.). 1990. Numerical Modeling of Seismic Wave Propagation. Soc. Expl. Geophys., Geophysics Reprint Series, no. 13. 525 pp.
• KRYLOV VI. 2006. Approximate Calculation of Integrals (Dover Books on Mathematics). New York, Dover Publications. 368 pp.
• MITTRA R & LEE SW. 1971. Analytical Techniques in the Theory of Guided Waves. Macmillan, New York. 320 pp.
• PETRASHEN GI. 1978. Basis of mathematical theory of elastic waves propagation. In: Questions of Dynamical Theory of Seismic Waves Propagation, vol. 18, Leningrad State University Publ., Leningrad, Russia, 31-53.
• PETRASHEN GI, MOLOTKOV LA & KRAUKLIS PV. 1982. Waves in Stratified Homogeneous Isotropic Media. Nauka Publ., Leningrad, Russia. 323 pp.
• PRIIMENKO V, AMARAL LH, MITROFANOV G, MISSÁGIA RM & BUENO AD. 2007. Problemas inversos dinâmicos para estimativa de parâmetros dos reservatórios baseados na sísmica multicomponente. Convęnio Petrobras/UENF-LENEP/Fundaçăo Bio-Rio 460.018.027-1. Instrumento Contratual: 2100.0012181.05-4, 140pp.
• THOMSON WT. 1950. Transmission of elastic waves through a stratified solid. Journal of Applied Physics, 21: 89-93.
• URSIN B. 1983. Review of elastic and electromagnetic waves propagation in horizontally layered media. Geophysics, 48: 1063-1081.
• ZOEPPRITZ K. 1919. Erdbebenwellen VIII B, On the reflection and penetration of seismic waves through unstable layers. Göttinger Nachr., 66-84.

Datas de Publicação

• Publicação nesta coleção
  01 Dez 2009
• Data do Fascículo
  Jun 2009

Histórico

• Aceito
  23 Jun 2009
• Revisado
  18 Jun 2008
• Recebido
  18 Jun 2008