Modelagem e migração em profundidade 2D em meios com simetria polar local

Modelagem e migração em profundidade 2D em meios com simetria polar local

Marco Antonio Cetale Santos^I; Djalma Manoel Soares Filho^II; Paulo Léo Manassi Osório^III

^IPontifícia Universidade Católica do Rio de Janeiro (PUC-RJ), Rua Nísia Floresta, 73, Andaraí, Rio de Janeiro, RJ, Brasil. Tels.: (21) 9135-5272 / 3527-1633 - E-mail: [e-mail]

^IIPETROBRAS, Rua Primeiros Sonhos, 190/303, Jardim Guanabara, Ilha do Governador, Rio de Janeiro, RJ, Brasil. Tels.: (21) 9505-7889 / 3865-6985 - E-mail: [e-mail]

^IIIPontifícia Universidade Católica do Rio de Janeiro (PUC-Rio)

RESUMO

Este trabalho propõe uma técnica do tipo rotação de fase (PSM - Phase-Shift Method), para migração em profundidade de dados sísmicos para meios localmente transversalmente isotrópicos (LTI), nos quais a direção do eixo de simetria varia continuamente ao longo das camadas. Para a modelagem em meios LTI, onde cada ponto da malha tem suas características definidas pelas velocidades de fase P e SV, parâmetros de Thomsen, densidade e inclinação do eixo de simetria, desenvolveu-se um método baseado na solução da equação elástica da onda por diferenças finitas. Na separação dos modos de onda qP e qSV dos sismogramas, implementou-se um algoritmo baseado na solução da equação de Christoffel. A migração para cada família de tiro comum é realizada somente por PSM. Nos resultados das migrações usando reflexões dos tipos qP-qP e qP-qSV, os horizontes foram localizados precisamente e verificou-se que o processo é estável em relação à variação do eixo de simetria. O método proposto é para aquisições sísmicas multicomponentes podendo ser aplicado em dados sísmicos marítimos convencionais, como também em dados provenientes de aquisições do tipo cabo de fundo e cabo vertical. Como o método proposto se baseia em algoritmos que utilizam PSM a sua implementação paralela pode ser altamente eficiente.

Palavras-chave: migração, modelagem, rotação de fase, anisotropia, meios localmente transversalmente isotrópicos, diferenças finitas.

ABSTRACT

This paper shows a technique based on the phase-shift method (PSM) to implement pre-stack depth migration on locally transverse isotropic media (LTI), with the symmetry axis direction varies continually along the layers. For seismic modeling, a generalization of the finite differences method for the solution of the elastic wave equation was used. With this procedure, it was possible to accommodate seismic modeling on LTI media defined by six parameters at each grid point, i.e. , density, P and S wave propagation velocities along the local symmetry axis, Thomsen parameters and the direction of the local symmetry axis itself. In order to separate from the seismograms the qP and qSV wavefields, an algorithm based on the Christoffel equation was implemented. The migration for each common shot gather is implemented solely by phase-shift based algorithms, which means that not only the depropagation of the registered wavefield, but also the generation of the time matrices involved in the imaging condition were obtained in this manner for each set of parameters at each depth level. The migration results using qP-qP and qP-qSV reflections show that the horizons were located precisely, and that the process is stable in relation to the symmetry axis variations. The proposed method is for multicomponent seismic acquisitions and might be applied to marine seismic data using streamers, or Ocean Bottom Cables or vertical cables. Since the proposed method uses phase-shift algorithms, its parallel implementation can be highly efficient.

Keywords: migration, modeling, phase-shift, anisotropy, locally transverse isotropic media, finite differences.

INTRODUÇÃO

Os algoritmos de imageamento e inversão dos dados sísmicos usados com sucesso na prospecção petrolífera, que utilizam as abordagens tradicionais, baseadas em meios isotrópicos e ondas refletidas compressionais, têm tido resultados satisfatórios na busca de grandes acumulações, contudo têm-se mostrado inadequados na prospecção de trapas estratigráficas e estruturais mais sutis (Tsvankin & Thomsen, 1994), pois necessita-se de imagens de maior definição e resolução. Para isso meios anisotrópicos e propagação da onda cisalhante devem ser consideradas.

Algoritmos de imageamento de dados sísmicos que exploram anisotropias do tipo polar com eixo fixo (horizontal ou vertical) são mais realísticos em situações em que os horizontes de interesse estão localizados abaixo de espessas camadas de folhelho ou de rochas intensamente fraturadas. Contudo, em áreas tectonicamente perturbadas, nas quais as camadas de folhelho ou fraturas encontram-se curvadas, a suposição de um meio transversalmente isotrópico com eixo fixo (vertical-TIV, horizontal-TIH e inclinado-TII), pode não ser realística.

Aprimoramentos não são vistos somente na qualidade da imagem, mas também num maior grau de consistência entre os dados sísmicos e os de poço, conseqüentemente levando a uma interpretação geológica mais confiável. Tal consistência foi observada por Uzcategui & Mujica (1995) que mostraram a correlação do aspecto geológico de uma seção migrada em profundidade com a sísmica de poço e dados de poço, especificamente em um campo venezuelano. Em outro exemplo, Rimmer (1990) mostrou que o processamento padrão subestimou a profundidade de um reservatório de petróleo sob uma camada de sal e posteriormente verificou-se que tal erro foi provocado pela presença de uma espessa camada de folhelho. Notou também que o erro foi menor nas regiões sob camadas mais finas. Existem também provas que sugerem a melhoria em imageamentos sub-basálticos como na margem (plataforma continental) do norte do Atlântico no Reino Unido (Druzhinin et al., 1999). Outros casos que mostram a influência da anisotropia foram apresentados por Meadows & Abriel (1994) numa análise no Golfo do México e por Banik (1984) na bacia do Mar do Norte. Observa-se que mesmo a anisotropia pouco pronunciada produz um efeito significativo.

Ball (1995) faz um estudo detalhado do imageamento em profundidade de falhas inclinadas no Zaire. O maior incentivo é a fraca definição das falhas e o mapeamento dos reservatórios nas zonas de falhas em uma sísmica 3D. Embora o problema seja ampliado por baixos coeficientes de reflexão e múltiplas, o maior fator é a anisotropia, que é confirmado na inferência desta anisotropia pela comparação das velocidades intervalares de empilhamento e das velocidades de poço. Como forma de garantir esteresultado, foi feita uma inspeção baseada em VSP (Vertical Seismic Profile), dados de poço e estudos laboratoriais dos testemunhos. A inspeção proveu a localização, origem, magnitude e natureza desta anisotropia, que neste caso foi provocada pela presença de carbonato fraturado. A migração pré-empilhamento em profundidade foi processada com o método de Kirchhoff usando a estimativa desta anisotropia. A consideração da anisotropia levou a uma imagem significativamente melhor. Outros resultados e observações com relação ao problema da anisotropia podem ser encontrados em MacBeth & Lynn (2000).

Como mencionado anteriormente, os trabalhos consideravam as estruturas anisotrópicas com eixo de simetria fixo sendo definidas como TIV, TIH e TII. Neste trabalho, supõe-se que a simetria polar é local (meios LTI) para as camadas de folhelho e fraturas, como sugerido em Thomsen (2002), isto é, cada ponto possui um eixo de simetria particular. Esta suposição permite modelar camadas anisotrópicas com variação do eixo de simetria ao longo da mesma. Uma possível estrutura que gere tal anisotropia seria uma camada de folhelho dobrada tectonicamente.

Neste trabalho, a parametrização das camadas anisotrópicas é feita para cada ponto do meio, sendo necessários os valores das velocidades de fase P e S na direção do eixo de simetria local, densidade, parâmetros de Thomsen (1986) e e d e a direção do eixo de simetria em relação a um sistema global de referência. Utilizando esta parametrização, desenvolveu-se um algoritmo para a geração de dados sintéticos considerando meios LTI. A implementação deste algoritmo é baseada na técnica de diferenças finitas (Zahradník & Hron, 1992), originalmente válido para meios elásticos isotrópicos. No algoritmo, usam-se aproximações de segunda ordem nas derivadas parciais, os parâmetros são introduzidos através de integrações ao longo das linhas da malha, utilizam-se bordas laterais e inferior de absorção (Cerjan et al., 1985) e aplica-se o formalismo do vácuo na borda superior do modelo (Zahradník & Priolo, 1994).

Nesta modelagem são gerados campos de deslocamento horizontal e vertical e para se processar a migração em profundidade, usando o método de rotação de fase, faz-se necessária a separação dos modos de onda. Desenvolveu-se um método no domínio da freqüência baseado em Dellinger (1991), para a separação das ondas qP e qSV em meios TI com eixo inclinado, quando a profundidade de referência na qual é realizada a separação encontra-se numa camada anisotrópica. Nos casos em que a separação dos modos de onda ocorre no interior de camadas isotrópicas, aplicam-se os operadores divergente e rotacional, como em Bulcão (2004), Sun & McMechan (2001) e Cunha Filho (1992).

METODOLOGIA

O processo para se obter os resultados sintéticos gerados neste trabalho está dividido na criação do modelo sintético, na modelagem, separação dos campos e na migração.

A etapa inicial consiste em definir um modelo 2D que acomode as características anisotrópicas das camadas. O modelo deve conter as velocidades de fase das ondas P (a ou v_p(0° )) e SV (b ou v_sv(0° )), a densidade r, os parâmetros de Thomsen (e e d) e o ângulo de inclinação f do eixo local em relação ao sistema global. Tais parâmetros são definidos em cada ponto da malha do modelo, isto é, o meio pode ser heterogêneo. Nascamadas isotrópicas do modelo, os parâmetros e, d e f são iguais a zero.

O algoritmo de modelagem elástica desenvolvido gera dados multicomponentes que se assemelham a dados reais a menos de ruídos e outras peculiaridades inerentes ao processo de aquisição. Os sismogramas das componentes elásticas devem ser separados nos campos de onda P e SV para se processar a migração por rotação de fase. Se a camada onde se adquirem os sinais for anisotrópica será necessário utilizar o algoritmo de separação aqui desenvolvido, que é muito eficiente em meios homogêneos.

Nos exemplos mostrados, a primeira camada é isotrópica, e neste caso, para a separação dos modos de onda P e SV, aplica-se o divergente e o rotacional, respectivamente. Cabe ressaltar que, caso o dado seja proveniente de modelagem acústica ou de aquisição marítima o sismograma gerado corresponde ao de onda P, não necessitando de separação de campos.

O próximo passo é efetuar uma filtragem FK e fazer o silenciamento da onda direta, porém, como o dado é simulado, estes processos foram substituídos pela retirada da onda direta a partir da subtração de um sismograma contendo somente a onda direta. Este sismograma é gerado por um modelo de dimensões iguais ao original e com velocidades de fase da primeira camada, sem conter interfaces. Assim o resultado desta subtração elimina o maior causador de ruídos na migração.

A migração usa técnicas de rotação de fase para resolver a equação acústica da onda, processando separadamente os sismogramas das ondas qP e qSV. Os algoritmos de depropagação da onda e o de cálculo dos tempos de trânsito são híbridos, isto é, calcula-se a rotação de fase em cada ponto da malha de formas distintas para os casos isotrópico e anisotrópico. Estes algoritmos usam os parâmetros do modelo a, b, r, e, d e f diferentemente da modelagem, pois eles processam separadamente as ondas qP e qSV.

Os exemplos deste trabalho utilizam como condição de imagem os tempos de trânsito calculados a partir da propagação das ondas P (ou qP) (Cetale Santos et al., 2005). Então ao se processar a depropagação dos sismogramas P e SV são obtidas imagens geradas com as ondas refletidas P-P e P-SV (ou qP-qP e qP-qSV). Também é possível obter a imagem das ondas refletidas SV-P e SV-SV (ou qSV-qP e qSV-SV), bastando para isso utilizar como condição de imagem os tempos de trânsito da onda SV (ou qSV) (Cetale Santos et al., 2005).

Geradas as imagens das ondas P e SV para cada tiro comum (PSDM) é necessário processar separadamente o empilhamento das mesmas. Para o empilhamento das ondas P faz-se a soma direta das imagens para obter a seção sísmica em profundidade final. No caso das ondas SV é necessário utilizar um processamento para corrigir a mudança de fase existente ao longo deum mesmo refletor, pois senão haverá cancelamentos mútuos no empilhamento. Este processamento não foi abordado neste trabalho.

MODELAGEM ELÁSTICA EM MEIOS LTI

O fenômeno ondulatório em meios elásticos anisotrópicos é regido pela seguinte equação:

ondet_ijdenota o tensor de esforços, f_i representa a i-ésima componente de densidade de forças externas, r a densidade, e u_i é a i-ésima componente do vetor deslocamento.

A relação entre o tensor de esforços e as derivadas espaciaisdo vetor deslocamento é dada pela lei de Hooke generalizada expressa por:

onde c_ijkl representa as componentes do tensor elástico que para um meio triclínico, que na notação Voigt (Carcione, 2001) é dado por:

Assim sendo, substituindo a Eq. (2) na Eq. (1), pode-se escrever a Eq. (1) em termos do vetor deslocamento:

Para o caso 2.5D, a Eq. (4) é simplificada suprimindo-se a componente u₂ e todas as dependências em relação ao segundo índice, isto é, explicitamente utiliza-se o seguinte sistema de equações, apresentado por Cetale Santos et al. (2003b).

As Eqs. (5a) e (5b) são válidas para meios 2.5D e fontes lineares infinitas perpendiculares ao plano formado pelos eixos horizontal e vertical, ou seja, x₁ e x₃.

Nos meios LTI, para cada ponto da malha, são definidos os valores das velocidades de fase das ondas P e S (a e b), os parâmetros de Thomsen e e d, a densidade r e o ângulo f, que define a direção do eixo de simetria polar no ponto. Com estes valores, excluindo a inclinação, obtêm-se as componentes do tensor elástico, em cada ponto, relativas ao sistema de coordenadas local, conforme as seguintes equações:

Desta forma, para se utilizar as Eqs. (5a) e (5b), expressas num sistema de coordenadas globais (com eixos horizontal e vertical, x e z respectivamente), realizam-se transformações de coordenadas em todos os pontos através da seguinte matriz de rotação (Thomsen, 2002).

em que f é o ângulo de inclinação do eixo de simetria em relação ao sistema global. As componentes do tensor elástico para um meio triclínico podem ser escritas como:

onde c¢_mnop(x_i) denota as componentes do tensor elástico para um meio transversamente isotrópico expressas no sistema local que explora sua simetria. O tensor c¢_mnop(x_i) na notação de Voigt é escrito como:

O tensor c_ijkl(x_i) da Eq. (8) representa um meio transversalmente isotrópico onde o eixo de simetria está inclinado em relação ao sistema global. Este tensor, quando escrito na notação de Voigt, apresenta 21 elementos diferentes. Porém, no caso 2.5D, somente 6 destes elementos são usados, não implicando em um acréscimo computacional significativo em relação a um meio que seja representado pela Eq. (9) que tem 5 elementos diferentes.

As derivadas parciais encontradas nas Eqs. (5a) e (5b) são discretizadas com aproximações de segunda ordem da forma apresentada em (Zahradník & Priolo, 1994). Com relação às derivadas temporais usou-se as aproximações centrais

onde U representa uma das componentes do vetor deslocamento, = U(i·Dx,j·Dz,n·Dt), Dx e Dz são os passos nas direções horizontal e vertical, respectivamente, e Dt é o intervalo de amostragem temporal. Em relação às derivadas parciais, com respeito às coordenadas x e z, temos aproximações distintas para as derivadas simples e mistas. Para as primeiras, ou seja, aquelas que envolvem uma única coordenada, x ou z, utilizam-se as aproximações:

onde a_S, a_N, a_E e a_W são integrações do parâmetro elástico realizadas ao longo da malha, respectivamente nas direções sul, norte, leste e oeste.

Para as derivadas mistas no interior do modelo, utilizam-se as aproximações ''curtas''

Na superfície de observação (j = 1), para atender ao formalismo do vácuo (Zahradník & Priolo, 1994), usa-se a seguinte aproximação ''longa''

onde a_SE, a_SW, a_NE e a_NW são integrações realizadas com orientação norte-sul para os pontos à sudeste, sudoeste, nordeste e noroeste do ponto do cálculo.

A outra derivada parcial mista relativa à Eq. (12b) é obtida analogamente à Eq. (13), porém definindo novos coeficientes a_ES, a_WS, a_EN e a_WN com orientação de integração leste-oeste (Zahradník & Hron, 1992).

No método proposto, os intervalos temporais e espaciais são definidos conforme Faria (1993), garantindo um processo estável e sem dispersão numérica.

Separação dos campos de deslocamento da onda em meios LTI

Atualmente, com a crescente necessidade da melhoria na resolução sísmica, necessita-se de processos de aquisição cada vez mais sofisticados. A aquisição multicomponente é uma destas tecnologias, sendo que para aproveitar estes dados torna-se necessário a separação dos modos de onda P, SV e SH. Uma das técnicas usadas na separação envolve os operadores divergente e rotacional para os meios isotrópicos (Clayton, 1981; Bulcão et al., 2003; Sun & Wang, 1999). No caso anisotrópico, pode-se utilizar a metodologia aqui estendida para meios LTI, que foi inicialmente descrita para meios transversalmente isotrópicos (TI) por Dellinger (1991).

Separação de campo isotrópico

Um campo vetorial arbitrário U(x,y,z) pode ser expresso por:

onde P é um escalar e S é um campo vetorial cujo divergente é igual a zero. Considerando duas dimensões ( x,0,z), S tem que ser representado por S = S1_y, onde S é um escalar e 1_y é o vetor unitário (0,1,0). Então ambos, P e S, podem ser tratados como escalares. Portanto, da Eq. (14), pode-se decompor um campo vetorial isotrópico U, nas suas componentes P e S, calculando-se Ñ·U = (ѲP) e ÑxU = (1_yѲS), respectivamente (Aki & Richards, 2002). Esta propriedade de separação de (Ñ·) e (Ñx) é explorada pelos algoritmos de migração e inversão elástica para o tratamento das ondas P e S (Bulcão et al., 2003).

Qualquer campo vetorial isotrópico ou não, pode ser separado como na Eq. (14), mas os potenciais resultantes representam ondas puras somente no caso isotrópico e homogêneo.

O cálculo do divergente (Ñ·) para obtenção do potencial P em duas dimensões é dado por:

Sua representação discretizada utilizando a aproximação de primeira ordem para a derivada é

onde Dx e Dz são os passos da integração espacial, u e w são as componentes horizontal (U_x) e vertical (U_z), respectivamente.

Usa-se o cálculo do rotacional (Ñx) para se obter o potencial S, Eq. (17), que por ser uma componente vertical, foi denominado de onda SV, dada por:

Sua representação discretizada é

onde Dx e Dz são os passos da integração espacial, u e w são as componentes horizontal (U_x) e vertical (U_z), respectivamente.

Outra forma de implementação dos operadores divergente e rotacional é no domínio da Freqüência. O operador divergente neste domínio é obtido pela transformada de Fourier da Eq. (15) e tem a seguinte representação:

onde () e () são as representações no domínio de Fourier de () e U, respectivamente, e

Reescrevendo a Eq. (19) de outra forma, temos a Eq. (20)

onde

é o número de onda em (k_x,k_z ) e é o vetor normalizado (k_x/k, k_z/k). Cada ponto (k_x ,k_z ) no domínio de Fourier representa ondas planas percorrendo na direção de , com freqüência espacial k. Em um dado ponto (k_x ,k_z ) há dois tipos de ondas, que são distinguidas pelas suas velocidades de propagação e pela direção da movimentação das partículas.No caso isotrópico, os tipos de onda são P e SV. Para a onda P as partículas se movimentam paralelas a e no caso SV perpendiculares a . Assim a Eq. (19) separa as ondas P, pois os movimentos das partículas são paralelos a , zerando as ondas SH que são perpendiculares.

O operador divergente no domínio de Fourier não é um separador puro, devido à presença do fator k. Este fator faz no domínio espacial a operação de derivação, amplificando as altas freqüências. Um operador que retire as ondas P e rejeite as SV sem alterar as freqüências pode ser construído pela transformada inversa de Fourier do operador (i ·), ao invés de utilizar o operador divergente,

= ik · . Entretanto, o divergente continua sendo muito utilizado, porque o operador i · apresenta uma descontinuidade em (k_x= 0, k_z= 0), onde é indefinido. Esta descontinuidade não é fatal, pois pode-se considerar um campo de onda que não tenha energia na freqüência zero, uma vez que os geofones são insensível a esta freqüência. Uma conseqüência desta descontinuidade é causar um operador no domínio espacial com o suporte menos compacto que o do divergente.

Caso anisotrópico

A principal mudança no caso isotrópico está na direção dos movimentos das partículas em relação à propagação dos modos de onda. Neste caso, os movimentos das partículas não são mais simplesmente paralelos e perpendiculares a . Agora esta direção de movimentação da partícula depende de uma expressão mais complexa dek_x e k_z.

Como no caso isotrópico, pode-se usar a característica da transformada de Fourier para se analisar separadamente a direção de movimentação das partículas para cada onda plana. Isso pode ser feito diretamente da equação elástica da onda no domínio de Fourier que é conhecida como equação de Christoffel (Carcione, 2001)

Simplificando a Eq. (20) para o caso 2D temos,

onde D é a matriz formada somente pelas freqüências espaciais normalizadas, l_x = k_x /k, l_z = k_z/k, C é a matriz das constantes elásticas que definem a propriedade do meio, r é a densidade do meio, w_n é a freqüência temporal para a solução n, k é o número de onda espacial e o_né a direção da movimentação da partícula para a onda plana (k_x,k_z) (Auld, 1973). Para um dado conjunto de parâmetros C, k e l, a Eq. (21) sempre tem duas soluções ortogonais o₁ e o₂ . Cada uma destas soluções corresponde a um modo de onda.

A separação dos modos de onda é feito por um algoritmo que determina a direção de movimentação das partículas no domínio da freqüência para cada modo de onda, e para um meio definido pela matriz de constantes elásticas C. Tendo-se a direção, pode-se facilmente separar o campo de ondas elásticas. Este algoritmo pode ser usado para encontrar o separador isotrópico e anisotrópico. Para o caso isotrópico usa-se diretamente o operador ik · na freqüência ou o divergente no espaço. O fluxograma deste algoritmo é mostrado na Figura 1.

Nas Figuras 2 e 3 são mostrados os campos de separação resultantes da aplicação do algoritmo da Figura 1 para o caso anisotrópico.

Para ilustrar a separação do campo, considera-se um meio anisotrópico homogêneo com parâmetros e = 0.3, d = 0.0, a = 3000 m/s, b = 1732 m/s, e densidade igual a 2460 kg/m³. Neste exemplo o modelo tem uma inclinação de eixo igual a 30° , e um valor extremo de e, que faz aparecer cúspides nas frentes de ondas propagadas (Fig. 4). Pode-se observar que a onda P na separação pelo método na freqüência está mais ressaltada (Fig. 6) do que na separação pelo divergente (Fig. 5).

PSDM por Rotação de Fase para meios LTI

Existe um grande número de trabalhos apresentando estratégias efetivas para o imageamento em profundidade em meios anisotrópicos, entretanto estes trabalhos levam em consideração meios homogêneos por parte ou eixo de simetria vertical, como por exemplo: Kitchenside (1991), Faria (1993) e Uzcategui (1994). Porém, em Thomsen (2002) é mostrado que para alguns casos, a suposição de se usar um eixo vertical de simetria ao longo do meio pode acarretar em erros significativos.

Nesta seção será apresentada a técnica de migração pré-empilhamento em profundidade (PSDM), que contempla meios transversalmente isotrópicos com qualquer inclinação do eixo de simetria, podendo ser definido de forma única em cada ponto do meio. Esta técnica está baseada no método de migração em profundidade por rotação de fase.

A equação de propagação da onda acústica tem como solução no domínio de Fourier

Após a aplicação da transformada inversa de Fourier obtém-se o campo p(x,z+Dz,t) na profundidade z+Dz, como mostrado por Gazdag (1978). Da relação de dispersão é obtido o módulo de k_z, e é aplicado o sinal negativo para uma solução ascendente, o que é dado por:

onde k_x, k_z , w e v são as freqüências espaciais, a freqüência temporal e a velocidade, respectivamente.

Para a modelagem, a solução da Eq. (22) propaga somente as ondas dos refletores para a superfície, usando um Dz com sinal negativo. No caso da migração Dz tem sinal positivo.

Esta metodologia é aplicada na migração pós-empilhamento, onde emprega-se a metade da velocidade de propagação para corrigir o tempo que a onda real leva para viajar da superfície até o refletor, somado ao tempo de retorno à superfície. A condição de imagem para esta situação é dada quando o tempo de depropagação for igual a zero.

Para a migração pré-empilhamento, faz-se a depropagação de cada sismograma separadamente e usa-se uma condição de imagem diferente da condição para a migração pós-empilhamento, dado pelo tempo de percurso da onda direta para cada ponto do meio. Assim a Eq. (22) toma a forma,

onde P é o campo que está sendo depropagado e T_D denota o tempo de percurso da onda direta.

A depropagação do campo para modelos LTI é processada através de um algoritmo híbrido, Cetale Santos (2003), ou seja, para os pontos isotrópicos da malha, k_z será computado de acordo com a Eq. (23), e para os pontos LTI da malha, k_z será estimado respeitando as seguintes equações:

onde q é o ângulo da normal à frente de onda e f é o ângulo de inclinação do eixo de simetria em relação ao sistema global.

O valor de k_zpara um meio isotrópico também poderia ser calculado por estas equações, mas o custo computacional para resolvê-las é maior do que o da Eq. (23).

Para melhor entender as Eqs. (25) e (26) parte-se da relação de dispersão:

que se mantém constante caso a velocidade não varie com o ângulo. Isso não ocorre nos meios transversalmente isotrópicos, pois a velocidade varia com o ângulo, tornando-se necessário usar as relações trigonométricas para se obter o valor de k_z. As relações das Eqs. (25) e (26) são retiradas do triângulo retângulo mostrado na Figura 7.

Para otimizar o processamento, ao invés de calcular as equações (25) e (26), estima-se a freqüência espacial k_z através da construção de uma tabela que relaciona sinq/ v(q-f) e cosq/ v(q-f) (Rousseau, 1997). Nesta construção varia-se o valor de q em intervalos pequenos de 0 a p rad. Esta tabela implicitamente relaciona os valores de k_x/w e k_z/w, segundo as Eqs. (25) e (26). Assim, tendo-se os valores de k_x e w é possível determinar k_z a partir da tabela.

Os valores de k_x e w que tornem o valor k_x/w fora dosvalores de sinq/ v(q-f) da tabela, representam o decaimento horizontal das ondas evanescentes, isto é, k_z é um número imaginário, e neste algoritmo quando ocorrem estes valores o campo é zerado. Isto significa, que não se levam em conta as ondas evanescentes.

As migrações para as ondas P-P e P-SV são computadas usando o algoritmo híbrido, onde a condição de imageamento para ambas é dada pelo tempo de percurso da onda de compressão.

As velocidades de fase são computadas pelas seguintes relações propostas por Thomsen (1986):

são as velocidades de fase, e e d são os parâmetros de Thomsen e q¢ = q-f, isto é, q' é o ângulo entre a normal da frente de onda (q) e o eixo local de simetria (f), como mostrado na Figura 8.

Na Figura 9 é mostrado o fluxograma do algoritmo de PSDM para meios LTI, utilizando a técnica descrida anteriormente. Neste algoritmo não são utilizadas interpolações de campos.

Para o cálculo do tempo de trânsito T_D é utilizado o mesmo esquema de rotação de fase híbrido (Cetale Santos et al., 2005). Entretanto, neste caso, se efetua uma modelagem descendente de uma fonte impulsiva localizada no mesmo lugar em que foi dado o tiro na modelagem LTI e usa-se o valor de k_z positivo.

Sendo o cálculo de T_D baseado na equação da onda em uma única direção, diminuem-se os erros que seriam provocados pelas reflexões, caso a equação da onda fosse completa.

Vale observar que na migração o ângulo q' é sempre o valor negativo do usado na modelagem. Na Migração a onda se propaga na direção de z crescente e tomando o refletor como referência, o ângulo inverte de sinal. Já para o cálculo de T_D não há esta necessidade, pois a referência está na superfície.

Esta técnica de migração LTI pode, também, efetuar migrações pós-empilhamento modificando-se a condição de imagem para T_D = 0 e usando-se nos modelos metade das velocidades de fase P e SV reais.

O algoritmo pode ser utilizado tanto para dados marítimos como para dados multicomponentes, pois trabalha com os modos de onda separadamente. Com algumas alterações, o mesmo pode ser usado em sistemas ortorrômbicos empregando as fórmulas de velocidade de fase apresentadas em Thomsen (2002).

Cabe ressaltar que com a técnica proposta não há neces-sidade de suavização dos parâmetros do modelo.

RESULTADOS

Nesta seção são apresentados três exemplos cobrindo aspectos do algoritmo quanto ao imageamento correto, a sensibilidade da migração, às variações do eixo de simetria e a sensibilidade da migração ao ruído aditivo.

Migração

O Exemplo 1 mostra o resultado da migração em profundidade de um meio LTI para uma família de tiro comum. Para tanto, gerou-se o modelo mostrado na Figura 10, onde é usada uma malha retangular de 400x800 pontos com intervalo espacial de 4m e seus parâmetros são definidos na Tabela 1. Este modelo tem uma estrutura complexa, que inclui uma falha e três camadas transversalmente isotrópicas com inclinações diferentes para os seus eixos de simetria. Para este caso, modela-se um tiro (família de tiro comum) com a fonte localizada no centro do modelo a 12m de profundidade, geofones cobrindo toda a extensão do modelo e intervalo temporal de simulação igual a 0.2 ms. Nas Figuras 11(a) e (b) são mostrados os sismogramas das componentes horizontal e vertical do campo de deslocamento, gerados na modelagem. Estes são processados para a retirada das ondas P e SV gerando os sismogramas das Figuras 11(c) e (d).

As Figuras 12(a) e (^b) mostram o resultado do PSDM anisotrópico, para ondas P e SV, respectivamente. Pode-se notar que as interfaces são recuperadas com precisão dentro da região bem iluminada pela geometria. No caso da imagem correspondente à onda SV, nota-se uma melhor resolução em relação ao da onda P, devido a seu menor comprimento de onda.

As Figuras 13(a) e (^b) mostram os resultados do PSDM isotrópico usando as velocidades de fase a e b, respectivamente. Note que no caso da imagem P, existe um pequeno erro de posicionamento em relação à interface horizontal na profundidade 1360 m, que não é perceptível no caso da imagem SV.

Sensibilidade do algoritmo de migração às variações dos parâmetros

No Exemplo 2 o modelo tem três camadas com interfaces planas como mostrado na Figura 14a. Neste caso a camada central é anisotrópica com inclinação do eixo de simetria igual a f = 20° e com e = 0.3 e o principal objetivo é recuperar a interface na profundidade de 1400 m. Os parâmetros completos deste modelo se encontram na Tabela 2. O PSDM foi processado para apenas um tiro, com uma fonte explosiva no centro do modelo e a uma profundidade de 12 m. Foi usada uma malha regular com 400x600 pontos e um intervalo de 4 m entre pontos consecutivos nas direções vertical e horizontal.

As Tabelas 3 e 4 mostram os erros de posicionamento da interface localizada a 1400 m de profundidade, relativos aos diferentes modelos de migração. Os erros foram obtidos pela diferença entre o posicionamento correto e o posicionamento do maior pico da forma de onda da imagem migrada para esta interface. Estes erros foram computados para 3 traços em 1100, 1200 (posição da fonte) e 1300 m. Os resultados da Tabela 3 foram obtidos com um processamento isotrópico, onde no primeiro teste é usada a velocidade vertical que é calculada a partir dos parâmetros anisotrópicos, e no segundo teste foram usadas as velocidades de fase a e b.

Na Tabela 4 temos os resultados obtidos nas simulações, com a mudança na inclinação do eixo de simetria, mantendo-se os valores de a, b, e, d e r usados na modelagem. Os resultados desta tabela demonstram a estabilidade em relação à variação do eixo de simetria.

Como mostrado em Cetale Santos et al. (2003a), para pequenas variações no ângulo do eixo de simetria, em torno de ± 5° para este exemplo, geram erros nas imagens da onda P pequenos e que se mantêm estáveis. Pode-se mostrar também que essas pequenas variações ocorrem para a onda SV. As variações nestes ângulos acarretam diferentes variações no posicionamento das interfaces para as ondas P e SV, o que era esperado, pois a velocidade de fase de cada modo de onda tem sua geometria própria.

Desta forma, os resultados obtidos supondo uma migração isotrópica são melhores do que os obtidos através da migração anisotrópica com erros na inclinação do eixo de simetria maiores que 5°.

Sensibilidade do algoritmo de migração ao ruído aditivo

O Exemplo 3 trata da sensibilidade do algoritmo de migração à introdução de ruído aleatório. Neste caso, são analisadas as seções migradas geradas a partir do sismograma de onda P,mostrado na Figura 11c, acrescido de diferentes níveis de ruído com distribuição normal. Os testes foram realizados pela adição de ruído branco, com relações sinal-ruído (RSR) de 10 dB, 1 dB e -10 dB, e pela adição de ruído colorido, com RSR de 1 dB,-10dB e -20dB a um sismograma sintético. O ruído colorido foi obtido pela filtragem do ruído branco por um filtro passa-baixa de resposta impulsional infinita, com freqüência de corte de 80 Hz. Assim, o ruído está na mesma faixa das freqüências do sinal sísmico.

A Figura 15(a) mostra o sismograma adicionado do ruído branco com RSR de -10 dB, e na Figura 15(b) um traço do sismograma no afastamento 800 m (representado pela linha branca mostrada em (a), e em (^c) o resultado da migração.

CONCLUSÕES

Foram efetuadas migrações separadamente das ondas P-P (qP-qP) e P-SV (qP-qSV) para meios LTI. Tais meios podem conter variações contínuas na inclinação do eixo. Para estas migrações, usou-se como condição de imagem os tempos de trânsito da onda P, mas nada impede que sejam usados os tempos de trânsito da onda SV (qSV), possibilitando assim as migrações das ondas SV-SV (e qSV-qSV) e SV-P (qSV-qP).

A técnica proposta, da forma que foi implementada, é mais precisa do que as do tipo PSPI (Phase Shift Plus Interpolation), pois não envolve interpolações. Isto é, por exemplo, se em uma determinada profundidade tiver 50 conjuntos de parâmetros, definindo 50 meios, a rotação de fase é feita para cada um destes meios e mesmo que todos os pontos da malha sejam diferentes, o algoritmo fará a rotação de fase para todos os pontos. Desta forma, o custo computacional aumenta de acordo com a complexidade do meio.

O algoritmo de separação dos campos de onda no domínio da freqüência é válido para meios anisotrópicos, porém está limitado a meios homogêneos. No caso de meios heterogêneos é possível traçar uma estratégia de solução que dividida o meio em partes homogêneas e use interpolação.

O algoritmo de modelagem desenvolvido permite observarque os meios anisotrópicos provocam grandes alterações nosdados sísmicos. Essas alterações incluem: mudança de fase até mesmo em refletores no sismograma de onda P; diferentes comportamentos para reflexões nas interfaces isotrópica-anisotrópica, anisotrópica-anisotrópica e anisotrópica-isotrópica; alterações de fase e mudanças de amplitudes do sinal em relação aos dados gerados em meios isotrópicos.

A maioria dos testes usou camadas com parâmetros de anisotropia fortes, em relação aos encontrados em bacias sedimentares (Thomsen,1986), e serviram para a verificação da precisão doalgoritmo. Nos testes apresentados mostrou-se que a anisotropia pode realmente alterar o posicionamento dos refletores, conforme previsto na literatura científica. Este fato é de suma importância para a indústria petrolífera, já que se um reservatório estiver sob camadas anisotrópicas pode haver variações significativas no seu dimensionamento.

A precisão do algoritmo de migração desenvolvido foi confirmada através dos resultados obtidos nos testes. Cabe ressaltar que com apenas uma família de tiro comum, usando um arranjo de receptores longo, foi possível recuperar praticamente todas as características do modelo. As imagens obtidas com os modos de onda SV (qSV) oferecem resolução superior às encontradas com as ondas P (qP) por apresentarem velocidades menores, isto é, uma fonte sísmica produz ondas SV com suporte espacial mais compacto que as ondas P.

A obtenção dos parâmetros de Thomsen medidos em laboratório deve levar em conta o sistema local de simetria polarconsiderando, assim, o acamamento ou fraturas existentes no testemunho. Em aplicações de dados reais os valores de parâmetros de Thomsen devem estar associados aos sistemas locais, pois a não observação deste fato pode comprometer o imageamento, como verificado na análise de estabilidade do processo. Concluiu-se que a estimativa incorreta do eixo de simetria de um meio anisotrópico, pode acarretar em resultado pior do que o obtido assumindo o meio como isotrópico.

AGRADECIMENTOS

Os autores gostariam de agradecer a PETROBRAS pelo suporte financeiro através do programa de pesquisa PRAVAP 19. Gostaríamos também de agradecer a Marcílio Castro de Matos e Felipe Prado Loureiro pelos seus comentários e sugestões.

Recebido em 8 dezembro, 2006 / Aceito em 28 junho, 2007

Received on December 8, 2006 / Accepted on June 28, 2007

NOTAS SOBRE OS AUTORES

Marco Antonio Cetale Santos. Graduou-se em Engenharia Eletrônica pela UERJ (1996). Obteve o seu título de Mestre em Processamento de Sinais na PUC-Rio (1999), tendo como tema da dissertação uma aplicação no sistema elétrico de potência. Em seguida, iniciou o doutorado e motivou-se com a geofísica onde aprendeu diversas técnicas de processamento de sinais inerentes. Obteve seu título de Doutor em Processamento de Sinais na PUC-Rio (2003). Hoje atua como Coordenador Técnico do LPS/PUC-Rio, trabalhando com o processamento sísmico, ajudando na orientação de alunos interessados no processamento de sinais aplicados a Geofísica e desenvolvendo projetos conjuntos com a Petrobras. Áreas de interesse: Modelagem e migração sísmicas, processamento de sinais geofísicos.

Djalma Manoel Soares Filho. Bacharelado e licenciado em física pela UFRJ (1977-1981), mestrado em física na área de teoria quântica de campos no Centro Brasileiro de Pesquisas Físicas (1982-1986). Obteve seu doutorado em geofísica, na área de tomografia sísmica, no convênio Petrobras - Universidade Federal da Bahia (1990-1994). Ingressou na Petrobras em janeiro de 1987, desenvolvendo trabalhos na área de processamento de dados sísmicos na unidade de negócios em Belém do Pará, no período de julho de 1987 a janeiro de 1990. Em junho de 1994, assumiu a coordenação de projetos em modelagem e migração sísmica em profundidade, no Centro de Pesquisas da Petrobras e atualmente ocupa o cargo de consultor sênior. Áreas de interesse: modelagem, migração, tomografia sísmica e ensino na área de geofísica.

Paulo Léo Manassi Osório (in memoriam ). Graduou-se em Engenharia Eletrônica pela Universidade Federal do Rio Grande do Sul (1967), tornou-se mestre em 1970 pela PUC-Rio e Ph.D. no ano de 1976 pela Universidade de Houston, ambos em processamento de sinais. Em 1976, iniciou suas tarefas didáticas na PUC-Rio. Em 1983, afastou-se para fazer o seu pós-doutorado na Universidade de Houston, onde desenvolveu pesquisas na área de Geofísica aplicada e atuou como professor associado. Foi professor da Escola Naval durante 25 anos. Como diretor do DEE, implantou o programa de doutorado em Engenharia Elétrica. Posteriormente readaptou o Laboratório de Processamento de Sinais (LPS/PUC-Rio) para projetos relacionados à Geofísica, dando início a trabalhos em parceria com a Petrobras. Durante 25 anos ministrou a disciplina Sinais e Sistemas para os novos geofísicos contratados pela Petrobras. Suas principais linhas de pesquisa eram análise sísmica, técnicas de processamento de sinais e Transformada Wavelet. Faleceu em 15 de agosto de 2006.

All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License

Av. Rio Branco, 156, sala 2510

20043-900 Rio de Janeiro RJ - Brazil

Tel. / Fax: [phone]

[e-mail] var _gaq = _gaq || []; _gaq.push(['_setAccount', 'UA-604844-1']); _gaq.push(['_trackPageview']); _gaq.push(['_setSampleRate', '70']); (function() { var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); })();

• AKI K & RICHARDS PG. 2002. Quantitative seismology. 2nd ed., University Science Books. 687 pp.
• AULD BA. 1973. Acoustic fields and waves in solids. 1: John Wiley and Sons, 435 pp.
• BALL G. 1995. Estimation of anisotropy and anisotropic 3-D prestack depth migration, offshore Zaire. Geophysics, 60: 1495-1513.
• BANIK NC. 1984. Velocity anisotropy of shales and depth estimation in the North Sea Basin. Geophysics, 49: 1411-1419.
• BULCÃO A. 2004. Migração reversa no tempo de dados sísmicosempregando operadores vetoriais. Tese de Doutorado, UFRJ, COPPE, Programa de Engenharia Civil, 349 pp.
• BULCÃO A, SOARES FILHO DM & MANSUR WJ. 2003. Migração reversa no tempo com operadores elásticos: imageamento com vários modos de ondas. Anais do 8° Congresso Internacional da Sociedade Brasileira de Geofísica, Rio de Janeiro, CD-ROM.
• CARCIONE J. 2001. Wave fields in real media: Wave propagation in anisotropic, anelastic and porous media. Pergamon Press, Handbook of Geophysical Exploration. Seismic Exploration, 31: 390 pp.
• CERJAN C, KOSLOFF D, KOSLOFF R & RESHEF M. 1985. A nonreflection boundary condition for discrete acoustic and elastic wave equation. Geophysics, 50: 705-708.
• CETALE SANTOS MA. 2003. Migração em profundidade por rotaçãode fase dos campos de onda qP e qSV em meios com simetria polarlocal. Tese de Doutorado, Departamento de Engenharia Elétrica, PUC-RJ, 97 pp.
• CETALE SANTOS MA, SOARES FILHO DM & OSÓRIO PLM. 2003a. Phase-Shift Depth Migration on Locally Transverse Isotropic Media: Stability in Relation to Symmetry Axis Variation at Each Medium Point. 73 Ann. Internat. Mtg., Soc. Expl. Geophys., Dallas, Texas, Expanded Abstracts, CD-ROM.
• CETALE SANTOS MA, SOARES FILHO DM, OSÓRIO PL & ROSA FILHO JC. 2003b. Modelagem sísmica em meios Localmente Transversalmente Isotrópicos. Anais do 8° Congresso Internacional da Sociedade Brasileira de Geofísica, Rio de Janeiro, CD-ROM.
• CETALE SANTOS MA, SOARES FILHO DM & OSÓRIO PLM. 2005. Cálculo da matriz de tempo de trânsito por rotação de fase em meiosLocalmente Tranversalmente Isotrópicos (LTI), Anais do 9° Congresso Internacional da Sociedade Brasileira de Geofísica, Salvador, CD-ROM.
• CLAYTON R. 1981. Wavefield inversion methods. SEP, 27: 81-93.
• CUNHA FILHO CA. 1992. Elastic Modeling in earth models. Ph.D. Thesis, Stanford University, 126 pp.
• DELLINGER JA. 1991. Anisotropic seismic wave propagation. Ph.D. Thesis, Stanford University, Department of Geophysics, 184 pp.
• DRUZHININ A, MacBETH C & HITCHEN K. 1999. Sub-basaltic depth imaging using stacking attributes. 69 Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, CD-ROM.
• FARIA EL. 1993. Modeling, migration and focusing analysis in transversely isotropic media, Ph.D. Thesis. The University of Texas at Austin, 196 pp.
• GAZDAG J. 1978. Wave equation migration with the phase shift method. Geophysics, 43: 1342-1351.
• KITCHENSIDE PW. 1991. Phase shift-based migration for transverseisotropy. 61 Ann. Internat. Mtg., Soc. Expl. Geophys, Dallas, Texas, Expanded Abstracts, CD-ROM.
• MacBETH C & LYNN BH. 2000. Applied Seismic Anisotropy: Theory, Background and Field Studies. Geophysics reprint series, No. 20,Series Editor Daniel A. Ebrom, 725 pp.
• MEADOWS MA & ABRIEL WL. 1994. 3D phase-shift migration in transversely isotropic media. 64 SEG meeting, Los Angeles, USA, Expanded Abstracts, 1331-1348.
• RIMMER DH. 1990. Model-based depth migration for subsalt structure in the Gulf of Mexico. 63 Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 881-883.
• ROUSSEAU JH Le. 1997. Depth migration in heterogeneous, transversely isotropic media with the phase-shift-plus-interpolation method. 67 Ann. Internat. Mtg., Soc. Expl. Geophys, Dallas, Texas, Expanded Abstracts, 1703-1706.
• SUN R & McMECHAN GA. 2001. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics, 66(5): 1519-1527.
• SUN R & WANG A. 1999. Scalar reverse-time depth migration of pre-stack elastic seismic data. 69 Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1457-1460.
• THOMSEN L. 2002. Understanding Seismic Anisotropy in Exploration and Exploitation, Distinguished Instructor Series, No. 5, SEG and EAGE, 238 pp.
• THOMSEN L. 1986. Weak elastic anisotropy. Geophysics, 51: 1954-1966.
• TSVANKIN I & THOMSEN L. 1994. Nonhyperbolic reflection moveout in anisotropic media. Geophysics, 59: 1290-1309.
• UZCATEGUI OJ & MUJICA DL. 1995. Anisotropic poststack depthmigration, eastern Venezuela. 65 Annual Internat. Mtg., Soc. Expl. Geophys. Expanded Abstracts, 1171-1174.
• UZCATEGUI O. 1994. Depth Migration in Transversely Isotropic Media with Explicit Extrapolators, Ph.D. thesis, Center for Wave Phenomena, Colorado School of Mines, Golden, Colorado 80401-1887. CWP-163, 140 pp.
• ZAHRADNÍK J & HRON F. 1992. Robust finite-difference scheme for elastic waves on coarse grids. Studia Geoph. et Geod., 36: 1-19.
• ZAHRADNÍK J & PRIOLO E. 1994. Heterogeneous formulations of electrodynamics equations and finite-difference schemes, Seismic waves in complex 3-D structures, Report 1, Dept. of Geoph., Charles University, 357 pp.

Datas de Publicação

• Publicação nesta coleção
  05 Dez 2007
• Data do Fascículo
  Jun 2007

Histórico

• Aceito
  28 Jun 2007
• Recebido
  08 Dez 2006