Uma implementação do método da refletividade para modelagem eficiente de dados sísmicos marítimos e de OBC

Uma implementação do método da refletividade para modelagem eficiente de dados sísmicos marítimos e de OBC

Sérgio Adriano Moura Oliveira

Laboratório de Engenharia de Exploração e Produção de Petróleo, Universidade Estadual do Norte Fluminense LENEP/UENF. Av. Brenand s/n, Imboassica, 27925-310 Macaé, RJ, Brasil. Tel.: [phone]; Fax: [phone] – E-mail: [e-mail]

RESUMO

Neste trabalho é demonstrado como o método da refletividade pode ser eficientemente implementado para simular aquisições sísmicas marítimas. Isto inclui o registro de pressão no caso da fonte e dos receptores (hidrofones) estarem localizados na camada de água, o deslocamento (ou velocidade de vibração) vertical e horizontal no caso dos receptores estarem plantados no leito do oceano, simulando aquisições multicomponentes tipo OBC (Ocean Bottom Cable) ou OBS (Ocean Bottom Sensor). Uma das vantagens do presente método é a facilidade de seleção dos eventos que farão parte da resposta, tais como múltiplas internas, múltiplas associadas à superfície livre e ondas convertidas no caso P-Sv. Os fenômenos da atenuação e da dispersão podem ser levados em conta facilmente considerando as velocidades complexas, já que o método é desenvolvido no domínio ω-p. O algoritmo proposto possui custo computacional bem menor quando comparado com métodos puramente numéricos como diferenças finitas ou elementos finitos. Sendo assim, no caso onde um modelo de camadas planas horizontais seja suficiente para descrever a subsuperfície, este torna-se uma ferramenta prática para tarefas como prever com bastante precisão e rigor a sensibilidade sísmica de um determinado reservatório com relação à substituição de fluidos ou prever como uma determinada mudança de propriedade petrofísica pode influenciar a resposta de AVO multicomponente.

Palavras-chave: Sísmica de reflexão, modelagem sísmica, método da refletividade, dados marítimos, OBC.

ABSTRACT

In this paper, a special implementation of the reflectivity method is proposed to simulate marine seismic acquisitions. The method was designed for efficient calculation of pressure field, in the case of source and receptors (hydrophones) located in the water column and the vertical or horizontal displacement (or velocity) in the case of simulation of multicomponent OBC (Ocean Bottom Cable) or OBS (Ocean Bottom Sensor) acquisition. One of the main advantages of the method is the easy choice of events that will be present in the response, like internal multiples, free surface multiples and converted waves. Moreover, modeling in the frequency domain makes it easy to handle absorption and dispersion in inelastic media. The proposed method has low computer cost when compared with pure numerical methods like finite difference or finite element. This makes it a valuable practical tool for accurate measurement of the seismic response under fluid substitution in a reservoir or to simulate the AVO response, when dealing with 1-D stratified media.

Keywords: Reflection seismology, seismic modeling, reflectivity method, marine data, OBC.

INTRODUÇÃO

Recentes avanços da sísmica, como aquisição multicomponente e sísmica 4D, vêm exigindo cada vez mais o suporte da modelagem sísmica. Esta é uma ferramenta fundamental para tarefas como prever a sensibilidade sísmica de um determinado reservatório com relação à substituição de fluidos ou prever como uma determinada mudança de propriedade petrofísica pode influenciar a resposta de AVO. Em muitas destas tarefas práticas, um modelo de camadas planas e horizontais já é o suficiente para elucidar várias questões ligadas à exploração ou à caracterização de reservatórios. Nestes casos, algoritmos de modelagem elástica baseada no método das diferenças finitas não são práticos, pois são extremamente caros computacionalmente e só se justificam no caso de modelos complexos de subsuperfície. Outra desvantagem desta técnica é a dificuldade de interpretação do sismograma sintético gerado, pois a solução numérica da equação da ondanão modela separadamente cada tipo de evento. O método do raio (e.g. Cerveny et al., 1977) é outra opção muito mais rápida computacionalmente, no entanto as implementações convencionais deste não levam em conta ondas refratadas e os efeitos da atenuação e dispersão.

No caso particular do modelo de camadas planas e horizontais o método da refletividade já é de muito conhecido, principalmente na sismologia ver, por exemplo, Fuchs & Müller (1971). Este também pode ser uma ótima opção para modelagem sísmica elástica. Este método é baseado na representação da transformada de Fourier do campo de onda em um meio estratificado através de integração ao longo do parâmetro do raio, ângulo ou número de onda horizontal.

Uma das vantagens do método da refletividade é a possibilidade de discriminar eventos, pois é possível modelar a resposta com ou sem determinado tipo de múltipla, onda de superfície, conversão, etc. Outra vantagem é que sendo um método no domínio da freqüência, os efeitos da atenuação e dispersão podem ser naturalmente acomodados considerando-se a velocidade complexa.

A implementação do método da refletividade depende docaso em estudo. Assim, na sismologia, os eventos de baixa freqüência relacionados às ondas de superfície são muito importantes. Já a sísmica de reflexão se concentra nas ondas de volume, assim, na integração de ondas planas, é necessário considerar somente os números de ondas reais. A seguir será visto como adaptar e otimizar o método da refletividade para simular as aquisições sísmicas onde os eventos de interesse são as reflexões P-P e P-Sv, estando a fonte localizada na camada de água e os receptores podendo estar na camada de água ou plantados no leito do oceano.

METODOLOGIA

O método da refletividade é uma solução formal para o problema da incidência de uma onda plana harmônica elástica em um pacote de camadas horizontais. A abordagem seguida neste texto é baseada em Müller (1985). Assim, as fórmulas para o cálculo do campo de onda de pressão na camada de água e do campo de deslocamento no fundo do mar são apresentadas usando-se o conceito de potencial escalar e vetorial. No método da refletividade a resposta do meio estratificado é expressa em termo dos coeficientes de reflexão generalizados que respondem pelas reflexões primárias, múltiplas associadas à presença da superfície livre, múltiplas internas e conversões sofridas pelas ondas no pacote de camadas. A resposta para fonte pontual é obtida através de integração ao longo do parâmetro do raio. Esta abordagem possibilita uma visão mais clara da física do problema e através de manipulações de suas equações é possível selecionar facilmente os tipos de evento de nosso interesse. Outra vantagem é que a estabilidade numérica do método é incondicional, o que nem sempre acontece em métodos que usam matrizes de propagação (ver, por exemplo, Kennett, 1983).

Neste trabalho serão apresentadas fórmulas que permitem calcular sismogramas sintéticos para as seguintes configurações de aquisição (ver Fig. 1): dado de pressão adquirido por hidrofones localizados na camada de água, que pode se referir tanto a um levantamento marítimo convencional com streamer ou a hidrofones localizados bem próximo do leito oceânico, simulando um dado de OBC-P e também dados das componentes de deslocamento vertical e horizontal de geofones localizados no leito oceânico.

CÁLCULO DO CAMPO DE PRESSÃO

Para calcular o campo tanto na camada de água como no fundo do mar serão usadas as funções potenciais escalar e vetorial. No caso em estudos, tais funções descrevem o campo devido a uma fonte pontual e isotrópica localizada na camada de água. Portanto a natureza do campo de onda em estudo é 3-D com simetria cilíndrica. O potencial escalar pode então ser expresso de acordo com a seguinte expressão (Aki & Richards, 2000):

Esta fórmula significa uma integração azimutal de ondas planas em coordenadas cilíndricas ao longo do parâmetro do raio p. Temos que a₀ = ( – p²)^1/2, sendo α₀ a velocidade de onda P na água, J₀ é a função de Bessel de ordem zero e F(ω) o espectro de freqüência do pulso sísmico. Como será visto em seguida, a função A₀ especifica a resposta do meio a uma onda plana harmônica com freqüência ω e parâmetro do raio p. O campo de pressão pode ser calculado a partir da Equação (1) através da relação pressão-velocidade em um fluido, que no domínio da freqüência é: ∇P = ω²ρ₀∇Φ₀, logo:

CÁLCULO DO CAMPO DE DESLOCAMENTO

Embora o que realmente seja medido em geofones de OBC seja a velocidade de vibração do fundo marinho, neste trabalho será considerado o deslocamento. Para que o movimento de vibração do fundo marinho seja calculado corretamente, serão considerados os potenciais dos seguintes eventos: onda P transmitida (Φ₁) e onda S convertida (Ψ₁) gerada na passagem da onda do fluido para o sólido, cujas fórmulas são respectivamente:

e também a onda P ascendente (Φ₂) gerada na interface logo abaixo do fundo do mar e sua posterior reflexão (Φ₂) e conversão (Ψ₂) geradas no fundo do mar (ver Fig. 1), cujas fórmulas são:

Nas Equações (3) a (7):

sendo α, β velocidades de onda P e S, respectivamente, J_k são as funções de Bessel de kº ordem; p é o parâmetro do raio; A_n e B_n são termos que incluem os efeitos de fase e amplitude devido a propagação, transmissão e reflexões das ondas da fonte até o receptor; e são os coeficientes de reflexão P-P e P-S, respectivamente, devido à incidência por baixo na interface água-sólido.

Para a onda S ascendente gerada por conversão P-S nainterface logo abaixo do fundo do mar (Ψ₃) e sua posterior reflexão (Ψ₄) e conversão (Φ₄) geradas no fundo do mar (ver Fig. 1), as fórmulas são:

onde, e são os coeficientes de reflexão S-S e S-P devido à incidência por baixo na interface sólido-água. Tais coeficientes dizem respeito a relações de amplitudes entre potenciais e,a rigor, devem ser obtidos solucionando o sistema de equações obtidas com condições de contorno específicas (ver, por exemplo, Berkhout, 1987).

Para introduzir os efeitos da atenuação e dispersão consideram-se as velocidades das ondas P e S complexas. A parte imaginária da velocidade é responsável pelo decaimento da amplitude enquanto a dependência com a freqüência da parte real é responsável pela dispersão. Ambas dependem do fator de qualidade do meio. As seguintes expressões para as velocidades podem ser derivadas do modelo do fator de qualidade constante (ver, por exemplo, Kjartansson, 1979)

onde α_n(ω) e β_n(ω) são as velocidades de fase das ondas P e S na n-ésima camada, para a freqüência angular ω. Os fatores de qualidade das ondas P e S são, respectivamente, q ^p e q ^s e ω₀ é uma freqüência angular de referência. Vale lembrar que a hipótese de que o fator q é constante na faixa de freqüências da sísmica de exploração tem sido confirmada em várias medições de campo (ver, por exemplo, McDonal et al., 1958).

O vetor deslocamento é calculado a partir dos potenciais através da seguinte expressão: = ∇Φ+ ∇×(0, Ψ, 0), que em coordenadas cilíndricas permite calcular as componentes vertical e radial:

sendo: Φ_T = ∑Φ_n e Ψ_T = ∑Ψ_n. Assim:

onde:

Por estas últimas equações fica claro que tanto a onda P como a S contribuem para as componentes vertical e horizontal de vibração.

Se não for necessário incluir o efeito da superfície livre, as funções A₀, A₁, A₂, B₁ e B₂ são dadas pelas seguintes expressões:

Tais funções incluem os termos de fase devido às trajetórias fonte – receptor e às refletividades e . Tais refletividades possuem o sentido de coeficiente de reflexão generalizado que responde pela resposta de onda plana de todo pacote de camadas abaixo da n-ésima interface. Uma fórmula iterativa, simples e estável para cálculo destes coeficientes é apresentada por Müller, 1985:

e

onde MT_i e MB_i são as matrizes de refletividade no topo e na base, respectivamente, da camada i; e são as matrizes de coeficientes de transmissão/reflexão, ascendentes e descendentes, da interface i + 1 e E é a matriz de deslocamento de fase. Esta fórmula iterativa é aplicada a partir da última interface e permite obter os coeficientes de reflexão generalizados na base ou topo da kª interface:

SELEÇÃO DE EVENTOS

As Equações (1), (2), (15) e (16) são a base da implementação do método da refletividade proposta neste trabalho. A seguir sãodesenvolvidas algumas fórmulas específicas para as refletividades e para as funções A_n e B_n de maneira a ser possível a seleção de eventos. Começando pelo caso onde as ondas convertidas não sejam de interesse. Para isto, considera-se as refletividades das conversões como sendo nulas: r_ps = r_sp = 0 e lembrando também que no caso considerado, r_ss = 0. Após algum desenvolvimento, partindo das equações matriciais (22) e (23),chega-se na seguinte equação escalar:

onde e são os de coeficientes de transmissão/reflexão PP, ascendentes e descendentes, da n-ésima interface. Pode-se selecionar apenas reflexões primárias ao expandir a Equação (25) em uma série de potências de r_pp e manter apenas o termo de primeira ordem. Assim a fórmula iterativapara r_pp fica simplesmente como:

onde a espessura da n-ésima camada é z_n. Se for importante incluir o efeito da superfície livre, deve-se considerar todas as reverberações possíveis na camada de água (Fig. 2). As reverberações incluem as múltiplas (Fig. 2a) e seus fantasmas (Fig. 2b) que possuem as seguintes expressões de onda plana:

O termo A₀ que inclui todos estes eventos possui então a seguinte expressão (sem contar a onda direta):

Para incluir o efeito da superfície livre no dado de OBC, repete-se o mesmo procedimento com relação ao somatório das ondas que se transmitem para o meio sólido toda vez que uma reverberação na camada de água incide no fundo do mar (ver Fig. 2):

Nestas equações z_s representa a distância da fonte ao fundo do mar, z_r é a distância dos receptores ao fundo do mar e z₀ é a espessura da camada de água.

OBTENÇÃO DE SISMOGRAMAS NO DOMÍNIO t-x

Para obtenção de sismogramas no domínio t-x é necessário realizar a integração ao longo do parâmetro do raio e da freqüência:

Nesta seção serão discutidos alguns pontos importantes desta etapa, bem como a estratégia adotada neste trabalho para resolver este problema. Para avaliar a Equação (30) é necessário antes de tudo escolher a ordem de integração. No chamado método espectral a integração ao longo do parâmetro do raio é realizada em primeiro lugar, obtendo-se o resultado intermediário no domínio ω-x. Isto pode ser feito através da transformada de Hankel, ou por integração numérica ao longo de um caminho conveniente no plano complexo. Em seguida, a integração ao longo da freqüência é realizada obtendo-se finalmente o sismograma no domínio t-x. Outra alternativa é realizar antes a integração ao longo da freqüência, obtendo o resultado intermediário no domínio τ-p e, após isto, integrar ao longo de p. Este método é conhecido na sismologia como método da vagarosidade (Fryer, 1980). A estratégia utilizada neste trabalho segue a linha do método espectral, ou seja, para cada freqüência, integrais do tipo a seguir são realizadas numericamente:

No entanto é bom que se entenda o comportamento do integrando no plano complexo para o caso particular em estudo. O termo

possui pontos de ramificação em , que também são pólos de onde partem linhas de ramificação ao longo de Re(p) Im(p) > 0 (ver Tsvankin, 1995). Os coeficientes de reflexão contidos na função F (p, ω) possuem termos a_n, que possuem pontos de ramificação em . No caso de incidência de um meio de menor para outro de maior velocidade, tais pontos são responsáveis pela presença de ondas refratas na solução. Caso tais coeficientes possuam pólos, estes serão responsáveis pela presença de ondas de superfície. No presente trabalho, é considerado que a velocidade seja complexa, de modo a levar em conta a atenuação, o que move estes pontos para o 1º e o 3º quadrante. Assim é possível fazer a integração ao longo do eixo real, em intervalos que englobem os parâmetros de raio de interesse, dada à posição de um determinado receptor. Para realizar esta integração é usado o método de Filon (Mallick & Frazer, 1987), pois este é mais preciso do que o método dos trapézios e permite a utilização de maiores intervalos de integração. Para evitar ruídos numéricos, é altamente recomendado que se evitem cortes abruptos de F (ω, p) nos extremos do intervalo de parâmetros do raio. Isto é conseguido simplesmente pela aplicação de uma função taper. A transformada rápida de Fourier (FFT) é usada para ir do domínio da freqüência para o tempo. Com relação a este procedimento, deve-se procurar se certificar que a energia do sinal irá decair significativamente até o final do tempo modelado. Caso contrário, eventos que deveriam aparecer em tempos mais longos aparecem precocemente. Deve-se lembrar que para modelar um sismograma cujo tempo máximo de registro seja T_max usando FFT, deve-se usar um intervalo de freqüência igual a 1/T_max. Para uma discussão mais detalhada acerca do método espectral ver, por exemplo, Kennett (2001).

EXEMPLOS

Para mostrar as possibilidades da metodologia proposta, elaborou-se um modelo contendo 40 camadas, cujos valores das velocidades, da densidade e dos fatores de qualidade são compatíveis com os encontrados na Bacia de Campos, em regiões de águas profundas. Este modelo engloba uma lâmina d'água de mil metros e mais quatro mil metros de sedimentos. Os valores de V_p, V_s e densidade na parte dos sedimentos do modelo podem ser visualizados na Figura 3. Os valores dos fatores de qualidade q ^p e q ^s variam respectivamente entre 50 e 90 e entre 30 e 60. Na água o valor de q ^p foi fixado em 300. As simulações para o caso do dado de streamer foram realizadas com os seguintes parâmetros: nº de receptores: 176, espaçamento entre receptores: 24 metros, profundidade da fonte: 10 metros, profundidade dos hidrofones: 10 metros, freqüência máxima do pulso: 50 Hz. Na modelagem dos sismogramas de OBC, foi usado o mesmo número de receptores e afastamentos do exemplo anterior. Adicionalmente foi simulado um sismograma OBC-P, que é um sismograma de pressão registrado a um metro do fundo do mar.

RESULTADOS

Na Figura 4 é exibido o sismograma simulando um dado de streamer computado levando-se em conta os efeitos da atenuação e dispersão e incluindo-se também as ondas convertidas. No sismograma da Figura 5 temos a resposta completa, ou seja, incluindo também as reverberações devido à presença da superfície livre. O sismograma OBC de velocidade de vibração vertical é exibido na Figura 6 e o de velocidade de vibração horizontal é exibido na Figura 7 (tais sismogramas foram obtidos simplesmente derivando-se com relação ao tempo os sismogramas de deslocamento). O sismograma OBC-P é exibido na Figura 8.

DISCUSSÃO

Ao incluir-se o efeito da atenuação, observa-se claramente uma perda de resolução temporal progressiva nos sismogramas calculados, o que é sempre observado em dados reais. É interessante notar também que as ondas convertidas estão presentes mesmo em dados de streamer. Estas são o resultado das conversões S-P de ondas ascendentes incidindo no fundo do mar (evento A no sismograma da Fig. 4). Tais eventos são mais visíveis para afastamentos maiores. Note que ao ativarmos a opção de superfície livre, as ondas convertidas parecem sermascaradas pelas reverberações (Evento B no sismograma da Fig. 5) que ocorrem na camada de água (comparar as Figs. 4 e 5). Com o objetivo de tornar as reflexões mais visíveis, foram aplicados ganho variável com o tempo e balanceamento intertraços nos sismogramas. Este procedimento acabou também por acentuar alguns ruídos numéricos devido à integração das ondasplanas, que aparecem discretamente nos resultados na forma de um evento horizontal no topo da primeira reflexão.

Quanto aos dados de OBC, observa-se claramente que asondas Sv se fazem notar no sismograma vertical e as ondascompressionais, por sua vez, se fazem notar no sismogramahorizontal (evento C na Fig. 7). Isto é esperado devido à conversão que ocorre na interface sólido-fluido. Observa-se que o sismograma OBC-P (Fig. 8) guarda uma semelhança muito grande com o sismograma OBC-Z, uma diferença marcante é que as mesmas múltiplas do fundo do mar e as múltiplas peg-leg, aparecem com polaridade invertidas nestes dois sismogramas. O que é esperado (ver, por exemplo, Brown et al., 2002).

CONCLUSÕES

A implementação do método da refletividade, segundo a estratégia proposta neste trabalho, permite otimizá-lo para modelar dado sísmico marítimo tanto de streamer, como de OBC. Note que, em termos de eficiência computacional, sempre que possível é melhor usar as equações iterativas escalares (25) e (26) para o cálculo dos coeficientes generalizados do que as equações matriciais (22) e (23). O tempo de computação dos sismogramas sintéticos, tendo em vista um número não muito grande de camadas, é pequeno o suficiente para que a modelagem possa ser usada de maneira iterativa e fazê-la competitiva frente ao método do raio. Nos exemplos expostos, o tempo de computação (considerando uma PC workstation padrão) foi de poucos minutos, enquanto o tempo para modelar o mesmo dado com os mesmos recursos computacionais usando-se diferenças finitas seria de algumas horas. No entanto, para um número muito elevado de camadas (maior que 500), um método puramente numérico para solução da equação da onda pode se tornar mais vantajoso. Como o presente algoritmo de modelagem é analítico (exceto pelas integrais para passar do domínio ω-p para o t-x) os resultados obtidos com este podem ser usados como padrão para aferir a precisão de esquemas de elementos/diferenças finitas. Por fim, vale ressaltar que a eficiência computacional do algoritmo proposto, sempre estará muito ligada a alguns detalhes de implementação. Isto depende, por exemplo, do cálculo numérico das integrais envolvendo as funções de Bessel. O código Fortran90 e o arquivo executável, denominado synsis, gerados a partir da metodologia apresentada, podem ser obtidos na seguinte página da internet: https://lenep.uenf.br/ ~ gir/software/software1.html.

AGRADECIMENTOS

O autor agradece as observações e os comentários feitos por dois revisores anônimos, que contribuíram para melhorar a qualidade deste trabalho.

Recebido em 28 fevereiro, 2008 / Aceito em 3 agosto, 2009

Received on February 28, 2008 / Accepted on August 3, 2009

NOTA SOBRE O AUTOR

Sérgio Adriano Moura Oliveira. Formou-se em Engenharia Elétrica pela Universidade Federal da Bahia (1993), doutor em Geofísica Aplicada pelo Programa de Geofísica Aplicada da Universidade Federal da Bahia – PPPG/UFBA (1998). Atualmente é Professor Associado do Laboratório de Engenharia de Exploração e Produção de Petróleo da Universidade Estadual do Norte Fluminense. É sócio da Sociedade Brasileira de Geofísica – SBGf e da Society of Exploration Geophysicists – SEG. Sua área de atuação é o método sísmico em geral, com especial interesse em modelagem da propagação de ondas sísmicas, inversão sísmica e análise de sinais.

• AKI K & RICHARDS PG. 2000. Quantitative seismology. University Science Books. 2 ed. 726 pp.
• BERKHOUT AJ. 1987. Applied Seismic Wave Theory. Elsevier. 371 pp.
• BROWN RJ, STEWART RR & LAWTON DC. 2002. A proposed polarity standard for multicomponent seismic data. Geophysics, 67: 10281037.
• CERVENY V, MOLOTKOV IA & PSENCIK I. 1977. Ray method in seismology. University of Prague Press. 212 pp.
• FRYER GJ. 1980. A slowness approach to the reflectivity method of seismogram synthesis. Geophys. J. Roy. Astron. Soc., 63: 747758.
• FUCHS K & MÜLLER G. 1971. Computation of synthetic seismograms with the reflectivity method and comparison with observations. Geophys. J. Roy. Astron. Soc., 23: 417433.
• KENNETT R. 1983. Seismic wave propagation in stratified media. Cambridge University Press. 342 pp.
• KENNETT R. 2001. The seismic wave field. Cambridge University Press. 366 pp.
• KJARTANSSON E. 1979. Constant Q Wave propagation and attenuation. J. Geophys. Res., 84: 47374748.
• MALLICK S & FRAZER LN. 1987. Practical aspects of reflectivity modeling. Geophysics, 52: 13551364.
• McDONAL FJ, ANGONA FA, MILLS RL, SENGBUSH RL, VAN NOSTRAND RG & WHITE JE. 1958. Attenuation of shear and compressional waves in Pierre shale. Geophysics, 23: 421439.
• MÜLLER G. 1985. The reflectivity method: a tutorial. J. Geophys., 58: 153174.
• TSVANKIN I. 1995. Seismic wavefields in layered isotropic media. Samizdat Press. 100 pp.

Datas de Publicação

• Publicação nesta coleção
  06 Abr 2010
• Data do Fascículo
  Set 2009

Histórico

• Recebido
  28 Fev 2008
• Aceito
  03 Ago 2009