Resumos

1. Introdução

O conceito de onda é um dos mais importantes e disseminados na Ciência. Tem papel de destaque em Física, Geofísica, Engenharias, Medicina e etc, sendo especialmente útil em diversas aplicações tecnológicas. Na Geofísica, em especial, este tema desempenha um papel central no Método Sísmico [¹1. P. Kearey, M. Brooks e Ian Hill, Geofísica de exploração (Oficina de Textos, São Paulo, 2009)., ²2. O. Yilmaz, Seismic Data Analysis. Processing, Inversion, and Interpretation of Seismica Data (Society of Exploration Geophysicists, Tulsa, 2000).].

O fenômeno físico da propagação de ondas mecânicas (ondas sísmicas) na interior da Terra é o fundamento maior da Sismologia, seja aplicada à exploração de hidrocarbonetos ou não. Na Sismologia Aplicada, são gerados artificialmente sismos de pequenas proporções e, assim como num terremoto, as ondas geradas atravessam diversos tipos de materiais líquidos e sólidos, como as rochas, mar e solos que compõe as camadas superficiais do nosso planeta, bem como nas suas camadas mais profundas, e.g., manto e núcleo. A energia que retorna às estações de medida pode ser utilizada para inferir informações valiosas das camadas superficiais da Terra ou mesmo da Terra como um todo, na sismologia aplicada ou na sismologia, respectivamente.

Este artigo é o segundo e último de uma série. Ele objetiva discutir conceitos físicos importantes – como a anisotropia e equações de propagação de ondas sísmicas em tais meios e a modelagem computacional – para leitores graduandos ou graduados em Física, especialmente professores do magistério superior, incentivando a utilização de métodos modernos dentro do escopo da Física e contribuindo para a disseminação de conceitos de modelagem numérica, assunto praticamente neglicenciado nos cursos de graduação em Física. Especificamente, visa apresentar os fundamentos teóricos e numéricos da propagação de ondas elásticas em meios anisotrópicos. Primeiro são deduzidas, de forma simples, as equações lineares da elastodinâmica para meios anisotrópicos gerais, utilizando-se um formalismo integral e depois tais equações são simplificadas para o chamado caso “pseudo-acústico” em meios com um grau de anisotropia simples (isotropia transversa vertical, do inglês, VTI, vertical transverse isotropy). Tais equações objetivam descrever apenas a parcela compressional da onda (onda P) em meios VTI. Por fim, discute-se a aplicação do algoritmo de propagação de ondas em meios anisotrópicos para o mapeamento de estrutura, mostra-se que neglicenciar os efeitos da anisotropia acarreta em problemas nas imagens obtidas, demonstrando-se a importância do correto tratamento da anisotropia em aplicações da equação da onda.

Como é bem disseminado, as ondas provenientes dos terremotos são as medidas principais de que se dispõe para o mapeamento do planeta Terra como um todos, na chamada Sismologia (Geofísica de Terra Sólida). No caso da Sismologia Aplicada (ou, simplesmente, Sísmica), similarmente as camadas de rochas sedimentares que compõem as bacias sedimentares podem ser mapeadas, em busca de estruturas que reúnam condições favoráveis para a acumulação de hidrocarbonetos. É importante observar que a sismologia, que tem como objeto de estudo a Terra como um todo, opera em frequências bem menores que Sísmica, sendo, portanto, capaz de “enxergar” estruturas de maiores dimensões do que a Sísmica, que, por sua vez, é capaz de delinear estruturas, na melhor das hipóteses, da ordem de algumas dezenas de metros.

Como mencionado, o Método Sísmico aplicado à exploração de hidrocarbonetos é empregado no interior de bacias sedimentares, o que pode ocorrer tanto em terra (onshore) quanto em mar (offshore). Ele é a principal ferramenta não somente para a identificação como também para o monitoramento de reservatórios. Pode-se dizer, sem exagero, que este método é o grande e maior responsável pela identificação dos reservatórios já encontrados, desde o início do século passado. Na verdade, está fazendo exatamente 100 anos que John Clearence Karcher provou que o método sísmico baseado em reflexões de ondas era capaz de identificar reservatórios de hidrocarbonetos [³3. B. Frehner, The Leading Edge 40, 404 (2021).]. Em tal método, como já mencionado, ondas sísmicas geradas artificialmente se propagam para o interior da região de interesse, onde sofrem reflexões e refrações diversas (descritas pela Lei de Snell), bem como difrações — quando encontram contrastes de impedância entre as camadas de rochas. Parte da desta energia, em especial aquela proveniente de reflexões, retorna à superfície da Terra, onde as ondas podem ser medidas por equipamentos próprios especialmente construídos.

O posterior processamento dos dados adquiridos em campo é realizado para inferir informações valiosas sobre a geologia das rochas sedimentares que compõe a bacia. Grande parte dos procedimentos realizados se refere ao processamento digital de sinais [²2. O. Yilmaz, Seismic Data Analysis. Processing, Inversion, and Interpretation of Seismica Data (Society of Exploration Geophysicists, Tulsa, 2000)., ⁴4. A.V. Oppenheim, R.W. Schafer e J.R. Buck, Discrete-time signal processing (Prentice hall Englewood Cliffs, Nova Jersey, 1989).] e.g., a atenuação de ruídos não desejados e soma de informações redundantes (empilhamento) com o objetivo de elevação da razão sinal-ruído.

Uma das últimas etapas do processamento é a chamada migração, onde objetiva-se obter um mapa dos refletores, i.e., das regiões associadas a descontinuidades de impedância do meio geológico [¹1. P. Kearey, M. Brooks e Ian Hill, Geofísica de exploração (Oficina de Textos, São Paulo, 2009)., ²2. O. Yilmaz, Seismic Data Analysis. Processing, Inversion, and Interpretation of Seismica Data (Society of Exploration Geophysicists, Tulsa, 2000).]. A migração é um problema inverso [⁵5. A. Tarantola, Inverse problem theory and methods for model parameter estimation (Society for Industrial and Applied Mathematics, Filadélfia, 2005).] simplificado no qual objetiva-se obter informações sobre o domínio de propagação, especificamente a posição dos refletores, a partir do conhecimento da solução da equação da onda em determinados pontos (em geral na superfície), informações estas medidas em campo, e das velocidades de propagação da onda, em geral de baixa frequência, as vezes chamado de macromodelo de velocidades.

A migração é uma inversão aproximada pois, diferente de um problema inverso estrito senso, nele é requerido como entrada o modelo de velocidades de propagação do meio. Na Geofísica Aplicada, tal modelo é obtido dentro do fluxo de processamento através de técnicas especiais denominada de análise de velocidade. Uma das técnicas mais modernas para análise de velocidades é baseada na resolução do problema inverso (não linear) linearizado tendo como base a equação acústica da onda, uma técnica muito dispendiosa computacionalmente e cara, enquanto técnicas mais antigas e baratas são baseadas em tomografia sísmica com os dados de superfície.

Assim como a inversão citada acima, as técnicas de migração mais poderosas disponíveis atualmente utilizam algoritmos de propagação de ondas para reposicionar o registro sísmico medido para as suas corretas posições em profundidade, uma vez que os efeitos decorrentes da propagação de ondas podem, de alguma forma, serem desfeitos lançando-se mão da equação da onda aplicada, por exemplo, no sentido inverso do tempo. Neste caso, o algoritmo de migração [⁶6. E. Baysal, D.D. Kosloff e J.W.C. Sherwood, Geophysics 48, 1514 (1983).] é denominado de Migração Reversa no Tempo (RTM, do inglês reverse-time migration), sendo exatamente o procedimento que inspirou este trabalho. Entretanto, a RTM foi proposta originalmente para o caso acústico e isotrópico, sendo utilizada apenas nas últimas décadas em caso práticos (da indústria). Como curiosidade, menciona-se que os algoritmos mais utilizados para migração sísmica são baseados em soluções de equações integrais da propagação de ondas em meios homogêneos (e.g., migração Kirchhoff), adotado em virtude da sua eficiência computacional.

Fisicamente, entretanto, o mapeamento geológico baseado em propagação de ondas é muito mais intuitivo e simples de compreender, embora demande mais recursos computacionais em virtude da necessidade de utilização de um algoritmos para a propagação de ondas e da consequente marcha no tempo que é necessária durante a aplicação da condição de imagem. Ele se fundamenta em dois princípios físicos básicos: o primeiro e mais conhecido é o que a equação da onda é invariante sob reversão temporal e o segundo é conhecido como princípio da reciprocidade. Tal princípio, afirma que as posições da fonte e do receptor podem ser trocadas sem que se altere o resultado da medida.

O método de mapeamento de estruturas baseado na propagação reversa de ondas, pode ser acústico (quando se consideram apenas ondas compressionais ou, simplesmente, ondas P) ou elástico (quanto as ondas cisalhantes ou transversais ou, simplesmente, ondas S também são utilizada). Neste sentido, a teoria mais simples capaz de explicar a propagação dos diferentes modos de onda existentes no registro sísmico (seja terremotos ou aquisições sísmicas com fonte artificial) é a chamada teoria da elastodinâmica. A princípio, esta teoria pode servir de base para construção de um algoritmo para mapeamento de refletores. Entretanto, este tipo de abordagem é mais complexa e requer passos adicionais de processamento, por exemplo, para a separação dos modos de onda P e S, complicação que o autor julga desnecessária para os objetivos deste artigo. Neste sentido, será discutido apenas a aplicação de formulações para meios VTI “pseudo-acústicos”.

De toda forma, na teoria elástica, o meio de propagação é suposto elástico, ou seja, uma vez que as tensões deixam de existir, cessam também as deformações e o meio retorna à configuração inicial, antes da aplicação da tensão. Os dois modos de onda (P e S), bem como as conversões entre eles e as ondas de superfície surgem das soluções destas equações.

A teoria acústica também é capaz de descrever o comportamento de ondas sísmicas, embora de forma menos completa, uma vez que tal teoria somente é capaz de levar em conta a propagação de ondas compressionais (ondas P). Na teoria do imageamento sísmico [²2. O. Yilmaz, Seismic Data Analysis. Processing, Inversion, and Interpretation of Seismica Data (Society of Exploration Geophysicists, Tulsa, 2000).], o modo de onda P é considerado o mais importante, sendo frequentemente utilizado isoladamente para imagear estruturas geológicas na exploração de hidrocarbonetos [⁷7. L.D. Bartolo, Rev. Bras. Ens. Fis. 43 (2021), unpublished.]. Para que este procedimento tenha sucesso, os dados são processados buscando-se realçar os eventos associados à propagação de ondas P, atenuando os eventos associados à propagação de ondas S. O primeiro motivo para utilizar prioritariamente as onda P na sísmica se deve às amplitudes das ondas compressionais serem maiores, frequentemente se sobressaindo no registro. Em segundo lugar, a consideração apenas de ondas compressionais simplifica o problema e reduz o número de parâmetros necessários para a migração, além de levar a algoritmos eficientes.

Em relação ao meio de propagação, independente da teoria utilizada, é sabido desde o século XIX que as rochas se comportam como meios anisotrópicos [⁸8. M.J.P. Rudzki, em: Über die Gestalt elastischer Wellen in Gesteinen, II. Studie aus der Theorie derErdbebenwellen (Bulletin International de l’Académie des Sciences de Cracovie, 1897).], ou seja, que suas propriedades elásticas variam com a direção. Ainda no início do século XX, a propagação de ondas elásticas em meios com simetria azimutal (transversamente isotrópicos ou simplesmente TI) foi considerada [⁹9. M.J.P. Rudzki, em: Sur la propagation d’une onde eelastique superficielle dans unmilieu transversalement isotrope (Bulletin International de l’Académie des Sciences de Cracovie, 1912).], embora apenas muito mais recentemente efeitos de anisotropia sejam levados em conta em algoritmos de processamento sísmico.

Neste trabalho, discute-se a propagação de ondas elásticas em meios anisotrópicos gerais e, em particular, em meios VTI. A apresentação é detalhada tanto no que se refere aos conceitos físicos quanto à construção dos algoritmos, que é apresentada de forma tutorial em material suplementar, através de pseudo-códigos. A partir de uma modelagem sintética de um meio complexo, de interesse geofísico, busca-se então verificar a influência de levar em conta a anisotropia ou não (considerando o meio isotrópico) na qualidade da imagem obtida.

A estrutura do artigo é a apresentada a seguir. Primeiro apresenta-se as equações da elastodinâmica linear para meios anisotrópicos e heterogêneos gerais. Em seguida, tais equações são simplificadas para descrever apenas o modo onda compressional em meios com simetria transversa vertical (VTI). Na sequência, soluções numéricas explícitas, utilizando-se o método das diferenças finitas (MDF), são desenvolvidas. Por fim, discute-se a aplicação do algoritmo de modelagem de ondas no mapeamento de estruturas geológicas em meios VTI, comparando-se o resultado com os resultados advindos da aplicação de um algoritmo isotrópico.

2. Equações gerais da elastodinâmica

Nesta seção, as equações que descrevem a propagação de ondas em meios elásticos heterogêneos e anisotrópicos gerais são apresentadas. Como será visto elas são consequência da aplicação direta de leis da física. Aqui será feita a obtenção das equações básicas de forma simples, tendo por base o importante livro texto de Wapenaar e Berkout [¹⁰10. C.P.A. Wapenaar e A.J. Berkhout, Elastic Wave Field Extrapolation: Redatuming of Single- and Multi-Component Seismic Data (Elsevier Science Ltd, New York, 1989).], onde a formulação integral da segunda lei de Newton é utilizada. Detalhamento teórico adicional sobre a elastodinâmica podem ser encontrados nas Refs [¹¹11. K.F. Graff, Wave Motion in Elastic Solids (Dover Publication, Inc, Mineola, 1991)., ¹²12. J.D. Achenbach, Wave Propagation in Elastic Solids. Applied Mathematics and Machanics (North-Holland Publishing Company, Amisterdam, 1973)., ¹³13. K. Aki e P.G. Richards, Quantitative seismology (University Science Books, Nova York, 2002)., ¹⁴14. R.S. Wu, R. Dmowska e V. Maupin, Advances in wave propagation in heterogeneous earth (Academic Pr, Cambridge, 2006)., ¹⁵15. I. Tsvankin, seismic signatures and analysis of reflection data in anisotropic media. Handbook of geophysical exploration (Elsevier, Amsterdã, 2001).], entre outros livros textos avançados de Física.

2.1. Conservação do momento linear

Considere uma região do espaço contendo um material sólido, cuja densidade é dada por $\rho (\overrightarrow{r},t)$, como mostrado na Fig. 1. Seja $\overrightarrow{v}(\overrightarrow{r},t)$ a velocidade de deslocamento das partículas, onde $\overrightarrow{r}$ é a posição e t o tempo. Considere, então, uma porção deste material de volume V, delimitada pela superfície S com normal $\widehat{n}$ como mostrado.

Figura 1
Região do espaço contendo o sólido de interesse.

A lei da conservação do momento linear estabelece que qualquer alteração do momento linear das partículas dentro do volume deve ser igual ao fluxo de momento linear entrando pela superfície, adicionado da resultante das forças que atuam nas partículas em V (2^a lei de Newton generalizada), ou seja, [¹⁰10. C.P.A. Wapenaar e A.J. Berkhout, Elastic Wave Field Extrapolation: Redatuming of Single- and Multi-Component Seismic Data (Elsevier Science Ltd, New York, 1989).]

sendo a força que atua no volume V dada por

A primeira parcela da Eq. (2) se refere às forças de contato que atuam na superfície da região considerada devido ao restante do meio (na região interna, pela Terceira Lei de Newton, estas forças se anulam), sendo ${\overrightarrow{\tau }}_n$ o vetor de força que atua em cada ponto da superfície (que depende da orientação da superfície). O vetor de força pode ser obtido por meio do tensor de tensão para uma direção qualquer. A segunda parcela representa a força externa aplicada em cada ponto do volume (força de volume), sendo $\overrightarrow{f}$ a densidade de força.

A única diferença da Eq. (1) para a expressão correspondente ao caso acústico [¹⁰10. C.P.A. Wapenaar e A.J. Berkhout, Elastic Wave Field Extrapolation: Redatuming of Single- and Multi-Component Seismic Data (Elsevier Science Ltd, New York, 1989).] é que no caso elástico está sendo utilizada uma tensão geral ${\overrightarrow{\tau }}_n$ que inclui tensões de cisalhamento enquanto no caso acústico a tensão é a mesma em todas as direções, sendo dada por $-p\widehat{n}$ (Princípio de Pascal).

Aplicando o teorema de Gauss na Eq. (1), obtém-se

com

Como o volume V é arbitrário, chega-se à

Esta equação vetorial expressa a conservação do momento linear (ou quantidade de movimento), sendo conhecida como equação não linear de movimento para meios elásticos. Considerando-se um sistema de coordenadas cartesianas, e utilizando a notação de Einstein, pode-se escrever facilmente a Eq. (5) como

onde τ_ij são as componentes do tensor de tensão.

Para diversas aplicações práticas, como é o caso da Geofísica, o interesse se restringe ao regime linear, sendo as equações lineares dadas por [¹⁰10. C.P.A. Wapenaar e A.J. Berkhout, Elastic Wave Field Extrapolation: Redatuming of Single- and Multi-Component Seismic Data (Elsevier Science Ltd, New York, 1989).]

2.2. Relações constitutivas

Para a descrição completa dos movimentos em um meio elástico, é necessário ainda levar em conta as chamadas relações constitutivas do meio, isto é, como se relacionam os tensores de tensão ${\tau }_{ij}(\overrightarrow{r},t)$ e de deformação $e_{ij}(\overrightarrow{r},t)$ (ambos tensores de segunda ordem). No regime linear, sua forma mais geral é dada pela lei de Hook generalizada, dada pela expressão

onde c_ijkl é o tensor de elasticidade que caracteriza completamente as propriedades elásticas do meio (a notação de Einstein é utilizada). A Equação (8) é a relação linear mais geral possível entre o tensor de tensão e o de deformação, sendo similar à $\overrightarrow{F}=k\overrightarrow{x}$. Neste sentido, considerando-se um prisma alongado na direção x, com comprimento L e área da seção transversa A, quando se traciona ele com uma força F na direção x, a expressão

onde Y é o módulo de Young (que, diferente da constante k da mola, depende apenas do material), descreve o comportamento do sólido no regime linear, i.e., a sua deformação Δx/L. No caso de sólidos, quando existem tensões mais complexas aplicadas (que não a tensão normal discutida no exemplo anterior), a princípio, o tensor de elasticidades tem 3⁴ = 81 componentes, por ser um tensor de quarta ordem. Entretanto, devido à existência de diversas simetrias física, este número se reduz substancialmente. Primeiro, devido ao equilíbrio de rotações, tem-se

uma vez que τ_ij = τ_ji. Além disso,

uma vez que e_ij = e_ij (ver Eq. 15). Por fim, assumindo deformações adiabática, pode-se demonstrar que [¹³13. K. Aki e P.G. Richards, Quantitative seismology (University Science Books, Nova York, 2002).]

As simetrias descritas acima, reduzem o número de coeficientes independentes do tensor de elasticidade para 21. Dependendo do grau de anisotropia, o número de coeficientes independentes se reduz ainda mais, sendo o caso mais simples aquele em que as propriedades elásticas independem da direção, isto é, o caso isotrópico. Neste caso, existem apenas duas componentes independentes que podem ser escritos em função do módulo de Young e da razão de Poisson. Na verdade, o caso acústico pode ser visto como o caso mais simples de elasticidade, onde existe apenas uma quantidade independente. Uma exposição elegante das diferentes formas do tensor de elasticidade, para diferentes graus de anisotropia, decorrentes de graus distintos de simetrias, é feita na Ref. [¹⁶16. L. Bos, P. Gibson, M. Kotchetov e M. Slawinski, Studia Geophysica et Geodaetica 48, 265 (2004).]. Neste trabalho, apenas será utilizado, sem demonstrar, o tensor de elasticidade para meios com simetria azimutal.

De qualquer forma, na presença de uma distribuição de fontes de tensão (ou deformação) externa, a Eq. (8) se torna

com

onde $h_{kl}(\overrightarrow{r},t)$ representa as componentes do tensor de deformação, referente à fonte externa aplicada. Substituindo a definição de deformação em termos de deslocamentos (definição deste tensor), i.e.,

na Eq. (13), chega-se à

Usando a simetria c_ijkl = c_ijlk (Eq. (11)), obtém-se

Por fim, considerando a velocidade $\overrightarrow{v}$ das partículas $\overrightarrow{v}={\partial }_t\overrightarrow{u}$, chega-se à relação tensão-velocidade linear

Tal equação, em conjunto com a Eq. (7), descreve completamente a propagação de ondas em meios elásticos anisotrópicos gerais. Trata-se de um sistema de equações hiperbólico de primeira ordem, em termos dos campos de velocidade e tensão, onde o meio é caracterizado tanto pela densidade como pelas componentes do tensor de elasticidade do meio, que podem variar ponto a ponto.

2.3. Equação em termos de tensões

Aqui, obtemos as equações de segunda ordem em termos apenas das tensões [¹⁷17. L.D. Bartolo, C. Dors e W.J. Mansur, Geophysical Prospecting 63, 1097 (2015).], que serão utilizadas na seção seguinte para chegar às equações pseudo-acústicas para o caso VTI.

Para chegar à equação desejada, divide-se a Eq. (7) por ρ (multiplica-se por b = 1/ρ), derivando-se em seguida em relação a x_l e multiplicando, por fim, por c_ijkl (contração de k e l). Com isso, segue

Derivando a Eq. (18) em relação ao tempo, tem-se

Somando, então, as Eq. (19) e (20), chega-se à

onde

3. Equações VTI pseudo-acústicas

As equações vistas anteriormente descrevem todos os modos de onda presentes em meios anisotrópicos quaisquer. Nesta seção, serão obtidas as equações que visam descrever apenas o modo de onda compressional (similar ao som), denominado de quase-P (ou simplesmente q-P) em meios VTI, i.e., meios anisotrópicos que possuem simetria azimutal vertical. As onda compressionais são denominada de quase-P, assim como feito para o modo cisalhante de onda quase-S porque não se caracterizam como perturbações nas direções longitudinal e transversal à direção de propagação, em virtude da diferença entre velocidade de fase e velocidade de grupo [¹¹11. K.F. Graff, Wave Motion in Elastic Solids (Dover Publication, Inc, Mineola, 1991).] no caso anisotrópico.

Será visto que tais equações para meios “pseudo-acústicos” VTI podem ser obtidas por meio de uma abordagem simples e direta [¹⁸18. L.D. Bartolo, L. Lopes e L.J.R. Lemos, Geophysics 82, T225 (2017).] tomando por base as equações vistas anteriormente e aplicando-se uma simplificação adequada, na qual o meio é considerado como “acústico”, embora anisotrópico. Empregá-se a palavra “pseudo-acústico” com o objetivo apenas de deixar claro que o interesse é em descrever o modo de onda compressional, muito embora fisicamente meios anisotrópicos não sejam acústicos, mas sim elásticos.

3.1. Notação de Voigt

Em virtude das simetrias existentes nos tensores de tensão, deformação e elasticidade (ver Eqs (10)–(12)), é possível reescrever estes campos com uma notação compacta. Isto é especialmente útil para o tensor de elasticidade, um tensor de alta ordem (quarta).

O tensor de tensão, por ser um tensor de segunda ordem, pode ser representado matricialmente por

onde, x = 1,y = 2,z = 3 e τ₂₃ = τ₃₂,τ₁₃ = τ₃₁,τ₁₂ = τ₂₁, totalizando 6 componentes independentes. Assim, uma forma conveniente de representar o tensor de tensão é a que segue:

onde, para tornar possível a representação utilizando um vetor, os índices i,j foram transformados em um índice n único, utilizando-se (notação de Voigt)

onde δ_ij é o delta de kronecker.

Utilizando a mesma notação — tendo em conta as simetrias nos índices originais ijkl vistas acima —, pode-se escrever o tensor de elasticidade como

Repare que as simetrias nos índices originais da tensão implicam que os novos índices são simétricos, evidenciando a existência de 21 propriedades independentes no caso anisotrópico geral, como mencionado acima.

3.2. Simetria VTI e parâmetros elásticos

Neste tópico, os parâmetros anisotrópicos que caracterizam a simetria vertical transversa (VTI) serão apresentados. Primeiro, entretanto, discute-se o que caracteriza este tipo de simetria, muito importante no caso da Geofísica, mas também em Física. Sua importância, no primeiro caso, reside no fato do depósito de sedimentos em bacias sedimentares se darem num sentido preferencial, ocorrendo um acamamento na direção horizontal na região mais central da bacial, como mostrado na parte (a) da Fig. 2. Com isso, nesta porção da bacia, gera-se uma rocha com simetria vertical transversa, i.e., um tipo de anisotropia na qual o módulo de elasticidade varia com o azimute mas é invariante por rotação em torno do eixo vertical de simetria. Na Figura 2, outro tipo de simetria azimutal é mostrado, na parte (b), na qual o acamamento de sedimentos ocorre na direçao vertical, sendo o eixo de simetria o eixo x: este meio é chamado de meio com simetria horizontal transversa (HTI, do inglês, horizontal transverse isotropy). Quando o eixo de simetria não é nem vertical nem horizontal, o meio é dito TTI (do inglês, tilted transverse isotropy).

Figura 2
Simetria VTI e HTI.

Agora, os parâmetros anisotrópicos que caracterizam os meios VTI serão apresentados. Aplicando-se as simetrias que caracterizam o meio (como a rotação em torno do eixo z), pode-se chegar a seguinte forma para o tensor de elasticidade em notação de Voigt [¹⁶16. L. Bos, P. Gibson, M. Kotchetov e M. Slawinski, Studia Geophysica et Geodaetica 48, 265 (2004).]:

Para chegar à aproximação acústica, naturalmente, deve-se exigir que c₅₅ = c₆₆ = 0 (parâmetros associados à velocidade de onda S), obtendo-se

Mais convenientemente, a utilização de parâmetros elásticos com significado físico (parâmetros de Thomsen e outros semelhantes) é muito útil. Os parâmetros de Thomsen (ϵ, δ e a velocidade c_Pz de propagação na direção vertical), de acordo com a definição dada na Ref. [¹⁹19. L. Thomsen, Geophysics 51, 1954 (1986).], podem ser escritos em termos de parâmetros elásticos, para o caso acústico, como

onde ϵ regula a diferença entre a velocidade da onda q-P nas direções vertical e horizontal: quando ϵ > 0, a velocidade é maior na direcção horizontal do que na direcção vertical; δ é responsável pela dependência angular da velocidade de fase das ondas q-P na vizinhança do eixo de simetria (vertical): quando δ > 0, a frente de onda se move mais rapidamente perto do eixo z em relação ao que ocorreria se o meio foi isotrópica. Este comportamento é exemplificado na Fig. 3. De toda forma, utilizando os parâmetros de Thomsen, pode-se escrever o tensor de elasticidade como

Figura 3
Parâmetros de anisotropia de Thomsen.

onde c_Pn (parâmetros que aparece no processamento sísmico) e c_PH são, respectivamente, a velocidade de sobre-tempo normal da onda q-P e a velocidade da onda q-P na direção horizontal, definidas como [¹⁵15. I. Tsvankin, seismic signatures and analysis of reflection data in anisotropic media. Handbook of geophysical exploration (Elsevier, Amsterdã, 2001).]

3.3. Equações em tensão e velocidade

Em notação de Voigt, a lei de Hook generalizada fica

No caso acústico VTI, é possível definir τ_H ≡ τ₁₁ = τ₂₂ e τ_V ≡ τ₃₃ (τ_ij = 0 se i ≠ j). Logo, tomando a derivada temporal da Eq. 34, pode-se escrever

Resta considerar as equações de movimento. A partir da Eq. 7 para o caso geral, chega-se facilmente à

As Equações 35-36 descrevem a propagação de ondas q-P em meios VTI 3D. As equações para meios 2D são semelhantes: a componente v_y desaparece, reduzindo-se de 5 para 4 equações. É importante notar que estas equações requerem que ε ≥ δ [²⁰20. T. Alkhalifah, Geophysics 65, 1239 (2000).].

3.4. Equações em termos de tensão

Eliminando as velocidades nas Eqs 35-36, obtém-se equações para o caso acústico VTI [¹⁸18. L.D. Bartolo, L. Lopes e L.J.R. Lemos, Geophysics 82, T225 (2017).]:

Uma vez que não se está interessados nas amplitudes das ondas para o mapeamento geológico, pode-se considerar a densidade constante ao longo do meio para reescrever as equações utilizando as velocidades:

É importante lembrar que a trajetória do raio associado à onda, que é o que de fato importa para o mapeamento de refletores, depende apenas das velocidades de propagação, de acordo com a lei de Snell. Além disso, é facil ver que as equações acima se simplificam para a equação acústica da onda no caso em que os parâmetros de Thomsen são negligenciáveis, pois

Logo as Eqs (39)-(40) se simplificam para a equação acústica da onda:

3.5. Fonte sísmica

De forma geral, uma fonte sísmica é qualquer aparato utilizado para gerar ondas sísmicas que irão se propagar através do meio. No formalismo desenvolvido neste trabalho, uma fonte geral é caracterizada pela densidade de força externa aplicada $\overrightarrow{f}$ e pela tensão externa aplicada g_ij. Por simplicidade, uma vez que rodamos apenas um exemplo onde a fonte é colocada no mar, considera-se apenas fontes explosivas, caracterizada por g_ij ≠ 0 apenas quando i = j e $\overrightarrow{f}=0$. Será adotado o chamado pulso de Ricker, que corresponde a uma g_ii Gaussiana 7, com g₁₁ = g₃₃ = i_V.

4. Solução numérica e aplicação ao “imageamento”

Nesta seção, será apresentada de forma direta uma formulação explícita simples a partir das Eqs 39–40 utilizando o método das diferenças finitas (MDF), bem como sua aplicação ao mapeamento e estrutura. Maiores esclarecimentos são encontrados no artigo que acompanha este [⁷7. L.D. Bartolo, Rev. Bras. Ens. Fis. 43 (2021), unpublished.]. Para conveniência do leitor, os algoritmo são discutidos em detalhes em material suplementar.

Antes de discutir a formulação numérica em si, lembramos que uma formulação explícita é tal que as soluções podem ser obtidas a partir do conhecimento das soluções em tempos anteriores. É importante mencionar também que o MDF é o método numérico mais simples capaz de fornecer soluções satisfatórias para equações diferenciais parciais. Ele divide o domínio do problema em pontos regularmente espaçados e utiliza aproximações (denominadas de diferenças finitas) para as derivadas, sendo um dos método numérico mais simples de programar. Ele fornece respostas muito eficientes em problemas de propagação de ondas, principalmente através da classe das chamadas formulações explícitas.

4.1. Formulação numérica adotada

As seguintes aproximações de diferenças finitas são utilizadas para obter a solução explícita do problema:

onde $F_{i,k}^n$ representa o campo f discretizado na malha de diferenças finitas, ou seja, em x = iΔx,z = kΔz,t = nΔt, sendo Δx e Δz os espaçamentos da malha e Δt o incremento temporal. As aproximações acima são, respectivamente, de segunda e quarta ordens.

As aproximações (e similares para a direção z) são substituídas nas Eqs 39–40 para o caso 2D, obtendo-se as equações de diferenças a partir da qual as soluções são geradas através da marcha no tempo explícita. Especificamente, denominando de (H,V) e (Z,X,N) as componentes discretas de (τ_H, τ_V) e (c_Pz,c_PH,c_Pn), isolando-se no primeiro membro o termo de tempo mais avançado, obtém-se

onde

Da mesma forma que a formulação acústica discutida na Ref. [²¹21. L.D. Bartolo, C. Dors e W.J. Mansur, Geophysics 77, T187 (2012).] e as formulações elásticas discutidas na Ref. [¹⁷17. L.D. Bartolo, C. Dors e W.J. Mansur, Geophysical Prospecting 63, 1097 (2015).], a formulação definida pelas equações acima é uma formulação explícita, como dito no início, uma vez que a resposta pode ser calculada explicitamente nos tempos futuros em função de tempos anteriores. Assim, especificamente, dadas condições iniciais conhecidas (em n = 0), a solução pode ser calculada de forma explicita no passo de tempo seguinte (em n = 1) em função de valores de $H_{i,k}^0$ e $V_{i,k}^0$ conhecidos no tempo nulo e de valores $H_{i,k}^{-1}$ e $V_{i,k}^{-1}$ no passo anterior (nulos pois anterior à aplicação da fonte).

De forma geral, conhecendo-se os valores dos campos H e V nos passos de tempo n e n−1, pode-se calcular todos os valores dos campos no passo de tempo n + 1 e este processo continua até que se atinja um determinado tempo final de análise desejado. Esquemas explícitos têm a vantagem de não necessitarem de solução de sistemas lineares a cada passo de tempo, um processo bem mais dispendioso computacionalmente. Os pontos de tensão normal discretos (sobre a malha de diferenças finitas) necessários para o cômputo de um ponto de tensão no tempo futuro definem o que se chama de estêncil do operador de diferenças finitas, mostrado na Fig. 4 para o caso em questão.

Figura 4
Estêncil do esquema acústico VTI (4E-2T).

Agora discute-se brevemente o “imageamento” de estruturas. Maiores detalhes devem ser procurados no artigo 1 [⁷7. L.D. Bartolo, Rev. Bras. Ens. Fis. 43 (2021), unpublished.], como mencionado anteriormente. Detalhes adicionais sobre a implementação computacional podem ser acessados no material suplementar. Como mencionado anteriormente, qualquer técnica de migração pressupõe o conhecimento do macromodelo de velocidades. Como se está trabalhando com dados sintéticos, o modelo de velocidade adotado é o modelo conhecido. Entretanto, para minimizar ruídos de baixa frequência inerentes ao método, utiliza-se este modelo suavizado, sendo realizada a suavização no inverso da velocidade (vagarosidade), como discutido em [⁷7. L.D. Bartolo, Rev. Bras. Ens. Fis. 43 (2021), unpublished.].

A primeira etapa da RTM consiste na propagação da onda no sentido direto do tempo, onde os tempos de chegada da onda direta em cada ponto do domínio são calculados e registrados na chamada Matriz de Tempo de Trânsito (MTT). Na segunda etapa, a energia do sismograma (obtido em campo ou sintético) que se deseja imagear é reintroduzida no domínio, através da aplicação dos valores medidos pelos receptores como fonte (na posição em que foram medidos) e propagando-se no sentido inverso do tempo, do final da análise até o inicio. Uma imagem para esta etapa é a seguinte: considerando-se a onda que propaga da fonte até os receptores com um filme, a segunda etapa nada mais é do que voltar o filme para trás, de forma que a onda vai retornando ao formato original.

Na segunda etapa, uma condição de imagem adequada deve ser aplicada. Ela é responsável por fazer com que a energia depropagada identifique o local em que ela foi gerada, seja ela uma reflexão ou difração. Em outras palavras, a condição de imagem é responsável por posicionar a energia depropagada na sua correta posição em profundidade, sobre os refletores. A condição de imagem estabelece que um ponto imagem P é todo aquele em que o tempo da onda depropagada é igual ao tempo da onda direta (dado pela MTT previamente calculada), sendo denominada de condição de imagem por tempo de excitação.

5. Exemplo numérico

No exemplo rodado, o mesmo do artigo 1, utilizou-se o modelo conhecido como Hess VTI (com a camada de água estendida), mas, neste caso, com a consideração da anisotropia VTI do modelo original (ver Figura 5). A geometria de aquisição (disposição das fontes e receptores para os diferentes tiros) adotada foi a mesma do artigo 1.

Figura 5
Modelo adaptador de Hess VTI. (a) velocidade de onda P na direção vertical, c_Pz, (b) ϵ e (c) δ.

O resultado final da migração, após o empilhamento dos 160 tiros, é mostrado na Fig 6 em comparação com a migração utilizando a formulação isotrópica [⁷7. L.D. Bartolo, Rev. Bras. Ens. Fis. 43 (2021), unpublished.]. Uma vez que os dados migrados são sintéticos, eles foram gerados utilizando a formulação “acústica” VTI, de forma a simular o comportamento VTI. Logicamente, ambas as imagens utilizaram estes mesmos dados. As imagens finais mostradas são apresentadas como saídas do algoritmo de empilhamento (no qual as images dos diferentes tiros foram simplesmente somada), não tendo sido aplicados quaisquer filtros.

Figura 6
Comparação da imagem final do modelo Hess após o empilhamento dos 165 tiros: (a) migração isotrópica e (b) migração VTI.

Pode ser visto claramente que a migração VTI produz resultados (Fig. 6-b) superiores e mais confiáveis que a migração isotrópica (Fig. 6-a), em especial, em uma área crítica: próxima às duas camadas delgadas que caracterizam o reservatório, na parte inferior direita do sal. Quando a anisotropia é tratada de forma adequada, a imagem ganha contornos mais nítidos e melhor foco na região. Além disso, pode-se ver também a atenuação de ruídos em comparação com a migração isotrópica. Enfatiza-se ainda que a imagem no caso isotrópico apresenta um erro de posicionamento de cerca de 80 metros no delineamento da base do sal.

6. Conclusões

Neste trabalho, um algoritmo explicito de modelagem acústica para meios verticalmente isotrópicos (VTI), baseado no método das diferenças finitas, foi aplicado para o caso de mapeamento de estruturas geológicas. O método de mapeamento de estruturas adotado foi inspirado no procedimento conhecido na literatura especializada como migração reversa no tempo pré-empilhamento. Trata-se de um dos métodos que melhores resultados tem fornecido, sendo utilizado especialmente para mapear regiões problemáticas para outros métodos, e.g., contornos de sal e baixo do sal.

As ideias fundamentais acerca da construção do algoritmo de modelagem em meios VTI foram discutidas, inclusive apresentado o algoritmo em detalhes em material suplementar, tendo sido a parte teórica abordada também com algum detalhamento, aquele que se julgou adequado aos objetivos do artigo. Foi discutido especialmente a questão do grau de anisotropia caracterizado pela simetria transversa. A aplicação da modelagem ao mapeamento foi discutida com detalhes. Viu-se que o método de mapeamento de estruturas em meios que apresentam comportamento anisotrópico pequeno (até um máximo de cerca de 0,2 para cada um dos parâmetros de Thomsen), embora significativo, demanda a utilização de equações de onda mais complexas, que levam em conta tais efeitos de anisotropia, para gerar resultados de melhor qualidade. Foi mostrado que estas equações de onda pra meios anisotrópicos recaem na equação acústica da onda no caso de meios isotrópicos, i.e., quando os parâmetros de anisotropia são negligenciáveis.

O resultado obtido no exemplo apresentado mostrou um erro de posicionamento de cerca de 80 metros no delineamento da base do sal quando se desconsidera a anisotropia do meio (mesmo esta sendo fraca), além de ruídos e mal posicionamento de refletores nos reservatórios. Tais problemas na imagem, no caso da aplicação a exploração de hidrocarbonetos, especialmente explorada neste artigo, podem fazer a diferença entre encontrar ou não petróleo em um poço exploratório, resultando no sucesso ou não de todo o projeto exploratório. Isto exemplifica bem a importância de se levar em conta efeitos de anisotropia em aplicações de propagação de ondas.

É importante mencionar que consideráveis desafios são introduzidos em aquisições reais, onde toda a informação (por exemplo, a velocidades de propagação e os parâmetros de anisotropia) deve ser extraída diretamente dos dados sísmicos, diferente da migração sintética aplicada neste trabalho, que utilizou os parâmetros conhecidos para a migração. Logicamente, a introdução de erros nestes valores leva a resultados tão piores quanto maiores os erros.

Por fim, salienta-se que os métodos e resultados aqui descritos podem ser extrapolados para casos de interesse em física, em especial casos que envolvem problemas inversos. Como discutido, o método desenvolvido é um problema inverso simplificado, onde informações de velocidades necessitam ser utilizadas. Frequentemente, quando aplicado em problemas reais de Geofísica, este tipo de mapeamento de estruturas é feito de forma iterativa, sendo o modelo de velocidades atualizado e o processo de migração refeito. Com isto, nos aproximamos ainda mais de um problema inverso. Entretanto, a grande vantagem do método descrito aqui é que o grau de complexidade matemática é consideravelmente menor do que em problemas inversos tradicionais, onde se busca minimizar uma função objetivo que quantifica a diferença entre dados modelados e medidos.

Desenvolvimento de algoritmos para propagação de ondas e mapeamento de estruturas em meios acústicos VTI

References

Datas de Publicação

Histórico