Propagação de ondas aplicadas ao mapeamento geológico: formulação acústica

Wave propagation applied to geologic imaging: acoustics

Resumos

O conceito de onda é recorrente nas ciências pura e aplicada, servindo de base para compreender desde terremotos, passando pelas telecomunicações e eletromagnetismo, até a difração de nêutrons, e.g., para citar um caso quântico. Assim, tal conceito desempenha um papel central em ramos tão diversos como Geofísica e Física Quântica, bem como em suas aplicações. A equação da onda é comumente estudada nos cursos de Física eminentemente utilizando soluções analíticas, não sendo possível obter respostas em domínios complexos e portanto estudar diversas aplicações. Este artigo objetiva apresentar detalhadamente, de forma simples e acessível, ricamente ilustrada, as ferramentas necessárias para obter soluções numéricas da equação da onda, dentro do escopo da Sismologia Aplicada, visando utilizá-las como base para um algoritmo de geração de imagens de estruturas geológicas em bacias sedimentares. Esta é uma bonita aplicação da Física, constituindo-se na base da exploração de hidrocarbonetos. A aplicação aqui descrita, conhecida na literatura especializada como Migração Reversa no Tempo, é aplicada a um modelo sintético, porém realístico. Os resultados demonstram a excelente qualidade da imagem obtida, sendo a técnica aplicada o estado da arte em termos de método de migração (mapeamento sísmico), atualmente muito em voga e utilizada pela indústria no pré-sal.

Palavras-chave
ondas acústicas; Métodos das Diferenças Finitas; Migração Reversa no Tempo; mapeamento de estruturas geológicas

The concept of wave is recurrent in pure and applied sciences, serving as a basis to understand and describe from earthquakes, through telecommunications and electromagnetic waves, to neutron diffraction, for instance, to cite a case involving Quantum Mechanics. Thus, this concept plays a central role in fields as diverse as Geophysics and Quantum Physics, as well as in their applications. On the other hand, the wave equation is commonly studied in Physics courses, eminently using analytical solutions, where it is not possible to obtain answers in complex domains and therefore study several applications, such as Applied Geophysics and Seismic Method. This article aims to present in detail, in a simple and accessible way, the necessary tools to obtain numerical solutions of the acoustic wave equation, within the scope of Applied Seismology, using them as the basis for an image generation algorithm (reflector mapping) of geological structures in sedimentary basins. This is a beautiful and elegant application of physics, forming the basis of hydrocarbon exploration nowadays, in which huge areas of the Earth’s surface are mapped with the objective of finding configurations in which there will be a potential commercial accumulation of oil. The application described here, known in the specialized literature as Reverse-Time Migration, is applied to a realistic synthetic model. The results demonstrate an excellent quality of the image obtained. In fact, the applied technique is state-of-the-art in terms of migration method (seismic mapping), currently very in vogue and used by industry, especially in pre-salt targets.

Keywords
acoustic waves; finite-diference method; reverse-time migration; geologic structure imaging

1. Introdução

O conceito de onda é recorrente em diversas áreas da Física e ciências afins, servindo de base para compreender e descrever fenômenos que vão desde a propagação de simples perturbações na superfície de um lago ou a propagação de terremotos através da Terra, até a propagação de luz, a difração de nêutrons entre outros fenômenos quânticos. Neste sentido, tal conceito desempenhou um papel central na nossa compreensão sobre a estrutura interna da Terra, de um lado, e sobre a estrutura da matéria, de outro. Desempenha ainda um papel importante em diversas aplicações científicas e tecnológicas no mundo contemporâneo.

Embora de importância inquestionável, a propagação de ondas é vista nos cursos de Bacharelado e Licenciatura em Física com uma abordagem estritamente analítica: apenas situações bastante simplificadas do ponto de vista prático são estudadas e perde-se a oportunidade de utilizar a equação da onda para uma abordagem introdutória de métodos numéricos, algo desejável na formação do cientista contemporâneo. Esta formação numérica é muito comum em áreas mais aplicadas como é o caso das engenharias e da Geofísica, embora infelizmente este não seja o caso da maioria dos cursos de Física pelo país; as cadeiras de física computacional ou similares raramente chegam a abordar métodos numéricos para resolução de equações diferenciais, especialmente EDPs.

Este trabalho visa justamente abordar o tema da resolução numérica da equação da onda, de forma introdutória mas cuidadosa – onde se discute a construção de algoritmos de modelagem e a aplicação deste ao mapeamento de estrutura utilizando ondas sísmicas – com o objetivo de contribuir para a disseminação de conceitos e técnicas numéricas entre professores.

Entre as inúmeras aplicações possíveis, este artigo escolheu abordar o mapeamento de estruturas utilizando a equação da onda, por ser uma aplicação muito atual e difundida em Geofísica, embora também apareça em Engenharia e Medicina e, potencialmente, em Física. Neste sentido, este trabalho tem como inspiração o método conhecido como Migração Reversa no Tempo (RTM, do inglês, reverse-time migration)[¹, ², ³], um dos métodos de mapeamento de estruturas geológicas mais poderosos utilizados em geofísica aplicada à exploração de hidrocarbonetos atualmente. Na RTM e em qualquer outro método de migração sísmica, como será discutido em detalhes, busca-se desfazer o efeito da propagação de ondas, reposicionando-se a energia sísmica proveniente de reflexões e difrações em interfaces onde existe contraste de impedância acústica.

A estrutura deste artigo é descrita a seguir. Na seção ², são discutidos conceitos importantes de sísmica de reflexão com o objetivo de tornar o texto o mais autoconsistente possível para físicos. Na seção ³, serão apresentadas as equações gerais que descrevem a propagação de ondas acústicas em meios materiais. Na seção ⁴, aborda-se a formulação numérica desenvolvida utilizando o método das diferenças finitas. Na seção ⁵, será abordada a aplicação do algoritmo descrito para o mapeamento de estruturas utilizando a propagação reversa de ondas, sendo o mesmo aplicado a um modelo sintético similar ao pré-sal brasileiro.

2. Conceitos fundamentais de Sísmica

Nesta seção, são discutidos alguns conceitos importantes sobre o Método Sísmico (ou, simplesmente, Sísmica) de forma bem direcionada aos nossos objetivos, ou seja, modelagem de propagação de ondas e sua a aplicação ao mapeamento de estruturas. Para uma abordagem completa e ainda introdutória, indicamos a Ref. [⁴]. Para uma abordagem mais detalhada e especializada, a Ref. [⁵] é a mais indicada.

O método sísmico de reflexão é o principal método da geofísica aplicada, no caso específico de exploração e caracterização de reservatórios de hidrocarbonetos (petróleo e gás natural), no sentido de que é o que recebe maiores investimentos por parte da indústria do petróleo. Por esta razão é também o método mais desenvolvido. Ele constitui-se de 3 etapas principais: aquisição, processamento e interpretação.

Na primeira, realizam-se in situ experimentos sísmicos cujo objetivo é fazer propagar ondas sísmicas na região de interesse, captando, em seguida, as ondas refletidas provenientes do meio ao longo do tempo. Na etapa de processamento, realizam-se procedimentos com o intuito de extrair as informações geológicas contidas nos dados coletados, gerando imagens da conformação geológica. De forma simplificada, pode-se dizer que o objetivo final desta etapa é obter uma imagem. Por fim, na terceira etapa, a imagem é interpretada (seja ela 2D ou 3D, como nas modernas aquisições e processamentos sísmicos), em conjunto com todas as informações geológicas e petrofísicas aferidas, visando identificar possíveis reservatórios.

As ondas que se propagam no interior de materiais sólidos em geral – como as camadas de rochas sedimentares que compõe as bacias sedimentares que abrigam os reservatórios de hidrocarbonetos – podem ser de dois tipos (Fig. 1): além das ondas longitudinais, acústicas, na sismologia chamadas de ondas principais (ou simplesmente ondas P), tratadas neste artigo, ocorrem também ondas cisalhantes (ondas secundárias ou ondas S), conversões entre estes tipos de onda e ondas de superfície. Claramente, as ondas S não se propagam em fluidos, uma vez que estes meios não apresentam resistência ao cisalhamento, como ocorre em material sólido.

Figura 1
Ondas de corpo: (a) onda P e (b) onda S.

Para descrever de forma mais detalhada a propagação de ondas em sólidos, é necessário lançar mão da chamada teoria da elasticidade (ou ainda de teorias mais sofisticadas), onde são corretamente tratadas tanto as ondas P quanto as ondas S, bem como as conversões de energia entre estes modos de onda quando da passagem das ondas através de interface onde há descontinuidade das propriedade físicas. Além disso, são corretamente descritas também as ondas de superfície, ondas que se propagam apenas próximo à superfície ou interfaces. As ondas P são mais rápidas que as S e ambos os modos de onda obedecem às leis usuais da refração e reflexão de ondas [⁶] (Lei de Snell):

(1) $$c_1\text{sen}{\theta }_1=c_2\text{sen}{\theta }_2$$

onde c₁ e c₂ são as velocidades de propagação da onda nos respectivos meios e θ₁ e θ₂ são os ângulos de incidência e transmissão (ou reflexão).

É importante mencionar que a teoria acústica, que será utilizada neste trabalho, descreve corretamente a propagação de ondas P na Terra, mesmo não sendo capaz de descrever as ondas S e as conversões. Assim, tal teoria, de fato, pode ser utilizada para o imageamento de estruturas geológicas desde que aplicada apenas à energia associada a ondas P medidas em campo, aplicação em que é empregada mais comumente na indústria, devido a maior eficiência computacional dos algoritmos acústicos.

Em relação aos tipos de levantamento realizados dentro do escopo do Método Sísmico, a chamada sísmica de reflexão é a mais comum. Ela objetiva coletar a energia de reflexões da onda sísmica que retornam à superfície após atingir regiões em subsuperfície com contraste de impedância acústica (I = ρc, sendo ρ e c, respectivamente a densidade e a velocidade de propagação da onda no meio). Outro tipo de aquisição é a chamada sísmica de refração, onde se objetiva captar a energia que sofre reflexão total (quando θ₂ = π/2, refração no jargão da área) em uma interface.

Em uma aquisição típica de sísmica de reflexão, uma fonte de ondas sísmicas é aplicada (usualmente próxima à superfície, em terra ou em mar), de forma que estas se propagam para o interior da área de interesse, ou seja, é gerado um sismo artificial de magnitude pequena. Quando encontra descontinuidades, a onda gerada pela fonte é refletida (como mencionado, proporcionalmente ao contraste de impedância acústica) e se dirige em direção à superfície, onde pode ser registrada e gravada para posteriormente ser processada e transformada em informação útil, especialmente para exploração de hidrocarbonetos. A fonte é detonada diversas vezes, varrendo, desta forma, imensas áreas e coletando grandes volumes de dados. A geometria de disposição das fontes e dos receptores influencia diretamente a qualidade final da imagem que poderá ser obtida. Outro aspecto fundamental é a redundância de informações, que possibilitará, no final das contas, aumentar a razão sinal-ruído e extrair informações sobre as velocidades de propagação do meio.

As aquisições feitas em mar aberto são operações realizadas de forma muito eficiente por navios próprios de aquisição sísmica. Os navios sísmicos atuais tem o tamanho de um campo de futebol e podem arrastar cabos de mais de 10 km. Nestas aquisições, são empregados arranjos de canhões de ar, dispostos próximos ao navio, como geradores de energia símica. Em navios de última geração, é possível utilizar 20 cabos ou um pouco mais (chamados de streamer), com distância entre eles de 50 m e comprimento de até 12 km. Nestes cabos estão dispostos os arranjos de hidrofones (receptores). Para se ter uma ideia de valores, se a distância entre tiros consecutivos em uma determinada aquisição é de 50 m, a velocidade do navio é de 5 nós (cerca de 2,5 m/s) e o tempo total de registro é de 6–10 s. Tipicamente, é possível varrer 100 km² ou mais em um turno de operação em águas livres. Uma aquisição de novos dados pode custar milhões de dólares.

As ondas registradas em superfície contém informações preciosas sobre a geologia em subsuperfície. Neste sentido, o conjunto de dados como um todo contém informações sobre as reflexões, refrações e difrações dos diversos modos de onda e de conversões destes modos de onda ao se propagarem na região de interesse. Mesmo em aquisições marítimas, onde são medidas as respostas apenas da componente compressional das ondas, há informações sobre as diversas conversões ocorridas em subsuperfície. Toda estas informações podem ser divididas em duas componentes, sendo a primeira relativa aos tempos de trânsito e a segunda relativa às amplitudes das ondas sísmicas. Os tempos de trânsito são utilizados para extrair informações chamadas de cinemáticas (no jargão geofísico), ou seja, informações relativas à posição dos refletores, as quais são influenciadas apenas pelo perfil de velocidades (de acordo com a Lei de Snell). Já as amplitudes podem ser utilizadas para obter informações de outra natureza (dinâmicas), relacionadas com os coeficientes de reflexão entre as camadas rochosas e consequentemente com diversas propriedades do meio, como a densidade do meio de propagação.

Em relação ao registro sísmico, após a geração artificial da perturbação pela fonte sísmica, o mesmo é feito por arranjos de equipamentos semelhantes a microfones, conhecidos como geofones (medido em terra) ou hidrofones (medido na água), dispostos, em geral, na superfície – no caso da chamada sísmica de superfície. Os geofones podem coletar dados multicomponentes enquanto os hidrofones coletam apenas informação da variação da pressão na água. Cada uma das medidas, coletada por um equipamento, é denominada de traço sísmico. Ao conjunto de dados obtidos dá-se o nome de sismogramas. Os sismogramas associados a cada tiro são denominados de painéis de tiro comum (common shot gathers). Tais registros contém informações da variação da pressão (no caso da água) ao longo do tempo para cada hidrofone que, por sua vez, são arranjados a diferentes distâncias da fonte (offset), como já mencionado.

Na segunda etapa, após a aquisição, os dados sísmicos são processados digitalmente com o objetivo final de se gerar uma imagem da região de interesse. Nesta etapa, são realizados procedimentos com o intuito de atenuar sinais não desejados e aumentar a razão sinal-ruído. Ao final, são realizados procedimentos objetivando fazer uma imagem dos refletores, colapsando as reflexões e difrações presentes nos sismogramas para suas corretas posições em profundidade (ou mesmo em tempo), processo chamado de migração sísmica.

O termo migração está associado às antigas técnicas geométricas (baseadas em traçado de raio, utilizando princípios como a Lei de Snell), baseadas em registros sísmicos feitos em papel. Buscava-se migrar as posições dos eventos observados nos sismogramas analógicos para suas corretas posições em tempo usando instrumentos típicos de desenho. Tais métodos de migração manual foram utilizados durante grande parte da primeira metade do séc. XX, sendo baseados em princípios geométricos simples, válidos apenas para configurações geológicas bem comportadas, como exemplificado para o caso de refletor plano (Fig. 2).

Figura 2
Migração de um refletor inclinado: em (a), é mostrada uma configuração geológica, formada por um trecho de refletor plano inclinado, sendo mostradas em f_i e r_i, respectivamente, as posições iniciais e finais das fontes e receptores e, em (b), mostra-se a “migração” dos eventos A e B para suas corretas posições em tempo (A e B).

3. A equação acústica da onda

Para desenvolver uma teoria que descreva o movimento longitudinal das partículas quando uma onda acústica se propaga no meio sob consideração, deve-se evocar os princípios fundamentais da conservação da massa e da conservação do momento linear, como por exemplo no texto de Wapenaar e Berkout [⁷]. Aqui, a dedução não será reproduzida, podendo o leitor se referir ao texto mencionado para a obtenção rigorosa ou mesmo a textos básicos de física, e.g., em [⁸]. As equações em termos do campo vetorial de velocidade de deslocamento de partículas $\overrightarrow{v}$ e do campo acústico de pressão p são dadas por:

(2) $${\displaystyle \frac{1}{\kappa (\overrightarrow{r})}}{\displaystyle \frac{\partial p(\overrightarrow{r},t)}{\partial t}}+\overrightarrow{\nabla }.\overrightarrow{v}(\overrightarrow{r},t)={\displaystyle \frac{\partial i_V(\overrightarrow{r},t)}{\partial t}}$$

(3) $$\rho (\overrightarrow{r}){\displaystyle \frac{\partial \overrightarrow{v}(\overrightarrow{r},t)}{\partial t}}+\overrightarrow{\nabla }p(\overrightarrow{r},t)=\overrightarrow{f}(\overrightarrow{r},t),$$

onde $\kappa (\overrightarrow{r})$ e $\rho (\overrightarrow{r})$ são, respectivamente o módulo de compressibilidade e a densidade do meio, $\overrightarrow{f}$ é uma força aplicada por unidade de volume e i_V é um parâmetro de injeção de massa utilizado em aplicações da Geofísica para gerar a fonte sísmica na água (meios offshore), descrito em [⁷] – ele representa a injeção de massa por unidade de volume do meio.

As equações (2)–(3) descrevem a propagação da onda acústica no regime linear, sendo o campo de onda completamente representado através das variáveis $p(\overrightarrow{r},t)$ e $\overrightarrow{v}(\overrightarrow{r},t)$, enquanto o meio de propagação é totalmente representado por $\rho (\overrightarrow{r})$ e $\kappa (\overrightarrow{r})$, sendo a velocidade de propagação da onda dada por

(4) $$c=\sqrt{\frac{\kappa }{\rho }}.$$

Pode-se manipular matematicamente as equação acima para obter uma que seja escrita apenas em termos do campo de pressão. Para tal, deriva-se temporalmente a Eq. 2 e divide-se por $\rho (\overrightarrow{r})$ a Eq. 3, tomando-se em seguida seu divergente, i.e.,

(5) $${\displaystyle \frac{1}{\kappa }}{\displaystyle \frac{{\partial }^{2}p}{\partial t^2}}+\overrightarrow{\nabla }.{\displaystyle \frac{\partial \overrightarrow{v}}{\partial t}}={\displaystyle \frac{{\partial }^2{i}_V}{\partial t^2}}\mathrm{e}$$

(6) $$\overrightarrow{\nabla }.{\displaystyle \frac{\partial \overrightarrow{v}}{\partial t}}+\overrightarrow{\nabla }.\left({\displaystyle \frac{1}{\rho }}\overrightarrow{\nabla }p\right)=\overrightarrow{\nabla }.\left({\displaystyle \frac{1}{\rho }}\overrightarrow{f}\right).$$

Subtraindo a segunda equação na primeira, chega-se à

$$\overrightarrow{\nabla }.\left({\displaystyle \frac{1}{\rho }}\overrightarrow{\nabla }p\right)-{\displaystyle \frac{1}{\kappa }}{\displaystyle \frac{{\partial }^{2}p}{\partial t^2}}=-F(t),$$

ou, multiplicando-se a equação acima por ρ,

(7) $$\rho \overrightarrow{\nabla }.\left({\displaystyle \frac{1}{\rho }}\overrightarrow{\nabla }p\right)-{\displaystyle \frac{1}{c^2}}{\displaystyle \frac{{\partial }^{2}p}{\partial t^2}}=-\rho F(t),$$

onde F representa a fonte sísmica, dada por

(8) $$F(t)=\frac{{\partial }^2{i}_V}{\partial t^2}-\overrightarrow{\nabla }.\left(\frac{1}{\rho }\overrightarrow{f}\right).$$

A equação (7) pode ser simplificada considerando-se meios em que o gradiente de densidades é desprezável. Primeiro, utilizando a identidade $\overrightarrow{\nabla }.(h\overrightarrow{\nabla }g)=\overrightarrow{\nabla }h.\overrightarrow{\nabla }g+h{\nabla }^{2}g$, pode-se reescrever a equação na seguinte forma

$${\nabla }^{2}p+\rho \overrightarrow{\nabla }\left({\displaystyle \frac{1}{\rho }}\right).\overrightarrow{\nabla }p-{\displaystyle \frac{1}{c^2}}{\displaystyle \frac{{\partial }^{2}p}{\partial t^2}}=-\rho F(t).$$

Logo, derivando o inverso da densidade, chega-se à

(9) $${\nabla }^{2}p-\frac{1}{c^2}\frac{{\partial }^{2}p}{\partial t^2}-\frac{1}{\rho }\overrightarrow{\nabla }\rho .\overrightarrow{\nabla }p=-\rho F(t).$$

Consequentemente, em regiões onde o gradiente de ρ é negligenciável, é válida a conhecida equação da onda

(10) $${\nabla }^{2}p-\frac{1}{c^2}\frac{{\partial }^{2}p}{\partial t^2}=-\rho F(t).$$

Em regiões onde o gradiente de densidade não seja desprezável, esta equação não é uma boa aproximação. Um caso típico onde isto ocorre é no fundo do mar, quando a onda passa da água para as primeiras camadas de sedimentos. Neste caso, embora, a princípio, a direção de propagação e os tempos de trânsito sejam bem descritos pela equação, o mesmo não se verifica com a amplitude. No jargão da área, a equação honra a cinemática da propagação, mas não a dinâmica.

3.1. 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 que se deseja estudar. Exemplos de possíveis fontes sísmicas para levantamentos terrestres em geofísica são: detonação de dinamites enterradas, queda de pesos e emprego de caminhões vibratórios.

No caso dos levantamentos marítimos, em geral, a fonte é um canhão de ar comprimido (airgun). Esta injeção de ar ocorre ao longo de um curto intervalo de tempo, tipicamente da ordem de um décimo de segundo (o que fornece frequências típicas de algumas dezenas de Hertz). Ela gera uma frente de onda esférica que se propaga pela água até atingir as primeiras camadas de rocha no fundo oceânico e continua a se propagar pelas rochas que compõe a subsuperfície. Ao atingir regiões com forte variação das propriedades físicas (contrastes de impedância), sofre reflexões, refrações e difrações, sendo uma parte desta energia enviada novamente para a superfície, onde podem ser medidas pelos hidrofones

A onda gerada por um airgun é de tipo P, isto é, puramente longitudinal, o único tipo de onda que se propaga em fluídos (como é o caso da água). Logicamente, isto não significa que inexistam ondas S em aquisições marítimas, uma vez que há conversão da onda P gerada pela fonte para onda S ao atingir as primeiras camadas rochosas em subsuperfície. Entretanto, neste trabalho apenas as ondas P serão modeladas, uma vez que a teoria acústica é utilizada.

A fonte deve ter certas características desejáveis: primeiro, é fundamental que tenham um espectro específico, onde se verifica a existência de uma frequência máxima presente, em segundo lugar, elas deve possuir uma frequência adequada para se verificar as feições na escala de interesse e que, por outro lado, permita que a onda se propague por distâncias da escala de cerca de uma dezena de quilômetros.

No formalismo utilizado neste trabalho, uma fonte geral é caracterizada pela distribuição $\overrightarrow{f}$ de densidade de força (forças de volumes) e pelo parâmetro i_V de injeção de massa. Como feito tradicionalmente na geofísica, será utilizada uma fonte dada pelo pulso de Ricker, o que corresponde a uma injeção de massa gaussiana [⁹], como será visto. O temo fonte utilizado na Eq. 10 é, desta forma, dado pela seguinte função

(11) $$F(t)=\left[2{\pi }^3{\left(f_c{t}_d\right)}^2-1\right]\exp \left[-{\pi }^3{(f_c{t}_d)}^2\right],$$

onde f_c é um parâmetro relacionado com a frequência de corte f_corte e t_d é o tempo defasado, utilizado para deslocar o início da aplicação da fonte, para que o máximo (ou mínimo, no caso) da função seja deslocado de zero para um tempo t₀ positivo de forma que a expressão seja praticamente zero no início da análise (em t = 0) e cresça suavemente (sem descontinuidade), conforme mostra a Fig. 3. As expressões para os parâmetros descritos acima são dadas por

Figura 3
(a) Fonte de Ricker, com f_corte = 60Hz, e (b) sua transformada de Fourier.

(12) $$f_c={\displaystyle \frac{f_{corte}}{3\sqrt{\pi }}}$$

(13) $$t_d=t-t_0$$

(14) $$t_0={\displaystyle \frac{2\sqrt{\pi }}{f_{corte}}}.$$

Como a Eq. 11 vai a zero apenas no infinito, na prática, utiliza-se tal expressão truncada no tempo, ou seja, a seguinte expressão é utilizada para o temo fonte:

(15) $$\begin{array}{cc}F(t)=\left[2{\pi }^3{\left(f_c{t}_d\right)}^2-1\right]\mathit{}\hfill & t\in [0,t_M]\hfill \\ \exp \left[-{\pi }^3{(f_c{t}_d)}^2\right],\hfill & \\ F(t)=0,\hfill & t\in (t_M,\mathrm{\infty }),\hfill \end{array}$$

onde

(16) $$t_M=2t_0={\displaystyle \frac{4\sqrt{\pi }}{f_{corte}}}.$$

Definindo o parâmetros i_V como (Gaussiana)

(17) $$i_V(t)={\displaystyle \frac{A}{2\pi {(\pi f_c)}^2}}\exp \left[-\pi {(\pi f_c{t}_d)}^2\right],$$

onde A é uma amplitude, introduzida por conveniência e, além disso, considerando a densidade de forças externas nula, ou seja,

(18) $$\overrightarrow{f}=\overrightarrow{0},$$

obtém-se o termo fonte visto acima. Com isto todos os termos fonte nas equações da acústica podem ser obtidos a partir das derivadas primeira e segunda de i_V, que por sua vez podem são dadas por

(19) $${\displaystyle \frac{\partial i_V}{\partial t}}=-A{t}_d\exp \left[-\pi {(\pi f_c{t}_d)}^2\right]$$

(20) $${\displaystyle \frac{{\partial }^2{i}_V}{\partial t^2}}=A\left[2\pi {\left(\pi f_c{t}_d\right)}^2-1\right]\exp \left[-\pi {(\pi f_c{t}_d)}^2\right].$$

Na Figura 4, são mostrados os gráficos do parâmetro de injeção de massa i_V e de sua derivada primeira para A = 1. Repare que este parâmetro é adimensional, por representar um volume de injeção de massa por volume do meio.

Figura 4
(a) Parâmetro i_V injeção de massa (função de Gauss) e (b) derivada primeira de i_V.

A transformada de Fourier da fonte é dada por

(21) $$s(f)={\displaystyle \frac{2f^2}{{\pi }^2{f}_c^3}}\exp \left[{\displaystyle \frac{-f^2}{\pi f_c^2}}\right],$$

o que fornece o conteúdo de frequências associados. Com isto, vê-se que a fonte é de banda limitada, sendo f_corte a máxima frequência presente (Fig. 3(b)).

O fato da fonte ser de banda limitada é importantíssimo pois, em primeiro lugar, as camadas rochosas por onde a onda se propaga funcionam como um filtro de altas frequências, dissipando energia nesta faixa. Em segundo lugar, para representar o sinal da fonte em termos discretos, existe uma frequência limite que pode ser representada dependendo do tempo de amostragem t_A utilizado na gravação do registro sísmico (como estabelecido pelo Teorema da Amostragem). Em terceiro lugar, o método das diferenças finitas sofre restrições em relação à máxima frequência que ele consegue aproximar sem que ocorra dispersão numérica, como será discutido posteriormente.

3.2. Condições iniciais e de contorno

Como é bem conhecido da teoria de equações diferenciais parciais, é necessário considerar-se condições de contorno e condições iniciais adequadas para que o problema tenha solução bem definida. Ou seja, trata-se de um problema de valor inicial e de valor de contorno.

Em relação às condições iniciais, nos problemas de Sísmica necessita-se apenas de condições triviais. Desta forma, as condições iniciais adequadas e utilizada aqui são de campos de pressão nulo, bem como de derivada temporal nula. Em relação às condições de contorno, classicamente dois tipos de condições de contorno são comumente utilizadas, nomeadamente as condições essenciais (de Dirichlet), com valor do campo prescrito, e as condições naturais (de Neumann), com valor da derivada do campo prescrita. Aqui utilizamos condições essenciais triviais (p = 0) em todas as bordas.

No contexto da Sísmica, onde o domínio seria melhor representado por um meio semi-infinito, aplicar p = 0 em todos os bordos do domínio numérico não é adequado, uma que que isto faz com que toda a onda incidente sobre as interfaces seja refletida. Apenas o topo do modelo deveria ter este comportamento, uma vez que se trata da superfície do mar nas simulações aqui realizadas. Nos demais bordos (nas laterais e na borda inferior), deve-se recorrer a estratégias especiais, chamadas de condições de contorno não-reflexivas, para que as ondas que as atinjam sejam atenuadas e não gerem artefatos numéricos em virtude do truncamento do domínio físico em um domínio computacional de dimensões adequadas.

Dito de outra forma, a necessidade de utilização de condições de contorno não-reflexivas se deve ao fato de os problemas que se deseja modelar serem na verdade problemas de domínios semi-infinitos, onde as bordas laterais e inferior devem ser simuladas como “transparentes”, evitando a degradação das soluções pelas reflexões nas bordas numéricas.

Uma estratégia simples e muito utilizada para atenuar as reflexões indesejadas nas bordas do modelo constitui-se na aplicação de camadas de absorção [¹⁰]. Esta zona de amortecimento é criada nos pontos próximos a borda do modelo, como mostrado na Fig. 5, sendo a função de amortecimento para o lado esquerdo, em função do ponto espacial n, dada por

Figura 5
A figura da esquerda mostra a exponencial aplicada em 100 pontos na região próxima à borda esquerda e a figura da direita representa a região de amortecimento.

(22) $$f_a(n)=Exp\{-{[a(N_P-n)]}^2\},$$

onde a é o fator de amortecimento, N_P é o número de pontos total da camada de amortecimento (em cada um dos lados) e n é o ponto espacial. No exemplo apresentado, utilizou-se a = 2,5×10⁻⁴ e N_P = 100. Como mostrado, em tais pontos o campo é multiplicado por um fator de decaimento exponencial ponto a ponto, o que reduz gradativamente a amplitude da onda ao se aproximar da borda do modelo. É importante chamar a atenção para o fato de que esta exponencial deve ser aplicada em pontos suficientes (da ordem de 80 a 100 pontos) e deve ser suave, como ilustrado. Outro métodos muito popular na área é denominado de PML (perfectly matched layers) [¹¹], método que modifica as equações na região próxima às bordas para que seja contemplada absorção de energia de forma perfeita (embora, após discretizadas as equações, a absorção não seja perfeita).

Outra estratégia que pode ser combinada com as camadas de amortecimento é a utilização de equações de onda de sentido único (one-way equations). Neste sentido pomos escrever a equação unidimensional da onda podemos fazer:

(23) $$\frac{{\partial }^{2}p}{\partial x^2}-\frac{1}{c^2}\frac{{\partial }^{2}p}{\partial t^2}=\underset{paradireita}{\underbrace{\left(\frac{\partial }{\partial x}+\frac{1}{c}\frac{\partial }{\partial t}\right)}}\underset{paraesquerda}{\underbrace{\left(\frac{\partial }{\partial x}-\frac{1}{c}\frac{\partial }{\partial t}\right)}}p=0,$$

onde identificamos cada um dos operadores diferenciais como responsável por fazer a onda pra pagar para cada um dos sentidos (direita ou esquerda), que dão origem a solução completa da equação da onda, a conhecida solução de D’Alembert. A partir da fatoração acima, podemos escrever as equações unidirecionais da onda:

(24) $${\displaystyle \frac{\partial p}{\partial x}}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{\partial p}{\partial t}}=0$$

(25) $${\displaystyle \frac{\partial p}{\partial x}}-{\displaystyle \frac{1}{c}}{\displaystyle \frac{\partial p}{\partial t}}=0.$$

A segunda equação é aplicada nos pontos da borda esquerda, o que tem o efeito de reduzir a reflexão da onda que atinge esta borda enquanto se propaga no meio 2D. A primeira equação é aplicada na borda direita e uma equação similar para a componente z é aplicada na borda inferior.

Para exemplificar a eficácia destas estratégias, é apresentada uma modelagem acústica para um meio de três camadas planas horizontais e homogêneas (Fig. 6). O efeito da aplicação isolada da equação one-way [¹²] é mostrada na Fig. 6(a) e o efeito combinado da aplicação desta e das camadas de Cerjan é mostrado na Fig. 6(b). Pode-se ver, que as reflexões nas bordas não são mais visíveis, uma vez que apresentam amplitude consideravelmente menor que as demais amplitudes presentes. Pode-se obter bons resultados aplicando-se apenas as camadas de Cerjan, como feito neste trabalho.

Figura 6
Condições de contorno simulando meios semi-infinitos, aplicando-se: (a) a equação one-way da onda isoladamente e (b) esta em conjunto com as camadas de amortecimento de Cerjan. Apenas na borda superior utilizou-se p = 0.

4. Soluções numéricas

Neste seção, primeiro serão apresentadas as ideias básicas do Método das Diferenças Finitas (MDF), com o objetivo de descrever em detalhes a formulação numérica explícita que será adotada. Para um aprofundamento do tema, reporta-se à vasta bibliografia especializada sobre o tema, em especial, ao excelente livro texto [¹³]. Outros livros sobre o MDF contém informação adicional e introdutória são [¹⁴, ¹⁵], podendo ser consultados outros bons livros textos mais gerais de métodos numéricos [¹⁶, ¹⁷, ¹⁸, ¹⁹].

4.1. O Método das Diferenças Finitas

O Método das Diferenças Finitas (MDF) é um dos métodos numéricos mais simples e intuitivos utilizados para resolução de equações diferenciais, podendo ser utilizado para obter soluções numéricas de uma vasta gama de equações.

O seu princípio básico é o de que as equações diferenciais são válidas ponto a ponto do domínio, podendo ser empregadas aproximações de diferenças finitas ponto a ponto desta malha para obter equações algébricas que irão fornecer soluções numéricas das equações. A Fig. 7 mostra a malha, com espaçamentos fixos Δx e Δz, respectivamente, nas direções x e z. No MDF, usualmente, também é adotado um incremento temporal fixo para a solução de equações dependentes do tempo. Desta forma, as funções e propriedades de interesse são consideradas apenas nestes pontos discretos P(x_i,z_j,t_n) = P(iΔx,kΔz,nΔt). Repare que o eixo z é apontado para baixo e será considerado o zero na superfície da Terra.

Figura 7
Malha simples de Diferenças Finitas.

São obtidas então expressões aproximadas para as derivadas nestes pontos discretos, sendo elas denominadas aproximações de diferenças finitas ou simplesmente diferenças finitas, podendo estas serem de diferentes ordens de precisão, dependendo do número de pontos utilizados, como será visto a seguir.

4.1.1. Operadores de diferenças finitas

Considere f(x) então uma função contínua, bem como suas derivadas. Pode-se desenvolver esta função em série de Taylor em torno do ponto x₀:

(26) $$f(x)={\sum _{n=0}^{\mathrm{\infty }}\frac{1}{n!}\frac{d^{n}f(x)}{d{x}^n}|}_{x=x_0}{(x-x_0)}^n.$$

Para obter as diferenças finitas que serão utilizadas nas formulações numéricas deste trabalho, o valor da função deve ser desenvolvido em série de Taylor em diferentes pontos x = x_i + lΔx, onde l ∈ Z, em torno de x_i, isto é,

(27) $$f_{i\pm l}=f(x_i\pm l\mathrm{\Delta }x)={\sum _{n=0}^{\mathrm{\infty }}\frac{1}{n!}\frac{d^{n}f(x)}{d{x}^n}|}_{x=x_i}{(\pm l\mathrm{\Delta }x)}^n.$$

Soma-se, então, as diferentes expressões obtidas a fim de cancelar os termos até a ordem de interesse da derivada em questão. Como será visto, quanto maior a ordem de aproximação, mais pontos precisam ser considerados. A seguir, serão obtidos os operadores de diferenças finitas em segunda e quarta ordens¹ de aproximação em Δx.

Assim, considere os seguintes desenvolvimentos (l = ±1):

(28) $$f_{i+1}=f_i+{\displaystyle \frac{\mathrm{\Delta }x}{1!}}{\displaystyle \frac{d{f}_i}{dx}}+{\displaystyle \frac{\mathrm{\Delta }{x}^2}{2!}}{\displaystyle \frac{d^2{f}_i}{d{x}^2}}+{\displaystyle \frac{\mathrm{\Delta }{x}^3}{3!}}{\displaystyle \frac{d^3{f}_i}{d{x}^3}}$$ $$\mathrm{}+{\displaystyle \frac{\mathrm{\Delta }{x}^4}{4!}}{\displaystyle \frac{d^4{f}_i}{d{x}^4}}+{\displaystyle \frac{\mathrm{\Delta }{x}^5}{5!}}{\displaystyle \frac{d^5{f}_i}{d{x}^5}}+\mathrm{\dots }$$

(29) $$f_{i-1}=f_i-{\displaystyle \frac{\mathrm{\Delta }x}{1!}}{\displaystyle \frac{d{f}_i}{dx}}+{\displaystyle \frac{\mathrm{\Delta }{x}^2}{2!}}{\displaystyle \frac{d^2{f}_i}{d{x}^2}}-{\displaystyle \frac{\mathrm{\Delta }{x}^3}{3!}}{\displaystyle \frac{d^3{f}_i}{d{x}^3}}$$ $$\mathrm{}+{\displaystyle \frac{\mathrm{\Delta }{x}^4}{4!}}{\displaystyle \frac{d^4{f}_i}{d{x}^4}}-{\displaystyle \frac{\mathrm{\Delta }{x}^5}{5!}}{\displaystyle \frac{d^5{f}_i}{d{x}^5}}+\mathrm{\dots },$$

onde df_i/dx é a derivada considerada no ponto x = x_i.

Fazendo-se a Eq. 29 menos a Eq. 28, obtém-se a expressão em segunda ordem de aproximação:

(30) $${\displaystyle \frac{d{f}_i}{dx}}={\displaystyle \frac{f_{i+1}-f_{i-1}}{2\mathrm{\Delta }x}}-{\displaystyle \frac{\mathrm{\Delta }{x}^2}{3!}}{\displaystyle \frac{d^3{f}_i}{d{x}^3}}-{\displaystyle \frac{\mathrm{\Delta }{x}^4}{5!}}{\displaystyle \frac{d^5{f}_i}{d{x}^5}}-\mathrm{\dots },$$

sendo o erro de truncamento da ordem de Δx², ou seja, é uma aproximação de segunda ordem.

Para obter aproximações de quarta ordem, é necessário considerar, além dos pontos x_i±1, os pontos x_i±2, utilizando procedimento similar. Primeiro, considera-se a Eq. 27 nos pontos l = ±2 e subtrai-se uma equação da outra, isolando a derivada, obtendo-se:

(31) $${\displaystyle \frac{d{f}_i}{dx}}={\displaystyle \frac{f_{i+2}-f_{i-2}}{4\mathrm{\Delta }x}}-{\displaystyle \frac{4\mathrm{\Delta }{x}^2}{3!}}{\displaystyle \frac{d^3{f}_i}{d{x}^3}}-{\displaystyle \frac{16\mathrm{\Delta }{x}^4}{5!}}{\displaystyle \frac{d^5{f}_i}{d{x}^5}}-\mathrm{\dots },$$

que inclusive pode ser obtida diretamente da Eq. 30 tomando os pontos i±2 (e fazendo-se Δx→2Δx). Para obter o operador com erro em Δx⁴, é necessário cancelar o termo em Δx² na expressão acima. Isto é obtido multiplicando a Eq. 30 por 4 e subtraindo a Eq. 31 dela, de onde segue

(32) $${\displaystyle \frac{d{f}_i}{dx}}={\displaystyle \frac{4}{3}}\left({\displaystyle \frac{f_{i+1}-f_{i-1}}{2\mathrm{\Delta }x}}\right)-{\displaystyle \frac{1}{3}}{\displaystyle \frac{f_{i+2}-f_{i-2}}{4\mathrm{\Delta }x}}$$ $$\mathrm{}+{\displaystyle \frac{4\mathrm{\Delta }{x}^4}{5!}}{\displaystyle \frac{d^5{f}_i}{d{x}^5}}-\mathrm{\dots }$$

ou, de maneira compacta,

(33) $${\displaystyle \frac{d{f}_i}{dx}}={\displaystyle \frac{-f_{i+2}+8f_{i+1}-8f_{i-1}+f_{i-2}}{12\mathrm{\Delta }x}}+O(\mathrm{\Delta }{x}^4).$$

Por outro lado, somando as Eq. 28 e Eq. 29 e isolando a derivada segunda, obtém-se

(34) $${\displaystyle \frac{d^2{f}_i}{d{x}^2}}={\displaystyle \frac{f_{i+1}-2f_i+f_{i-1}}{\mathrm{\Delta }{x}^2}}-{\displaystyle \frac{2\mathrm{\Delta }{x}^2}{4!}}{\displaystyle \frac{d^4{f}_i}{d{x}^4}}+\mathrm{\dots },$$

que é uma aproximação de segunda ordem para a derivada segunda. A aproximação de quarta ordem é obtida pelo mesmo processo visto para a derivada primeira, resultando em

(35) $${\displaystyle \frac{d^2{f}_i}{d{x}^2}}={\displaystyle \frac{4}{3}}{\displaystyle \frac{f_{i+1}-2f_i+f_{i-1}}{\mathrm{\Delta }{x}^2}}-{\displaystyle \frac{1}{3}}{\displaystyle \frac{f_{i+2}-2f_i+f_{i-2}}{4\mathrm{\Delta }{x}^2}}$$ $$={\displaystyle \frac{-f_{i-2}+16f_{i-1}-30f_i+16f_{i+1}-f_{i+2}}{12\mathrm{\Delta }{x}^2}}.$$

4.1.2. Propriedades das formulações de MDF

O MDF permite resolver equações diferenciais parciais (EDP), como é o caso da equação da onda, através da aplicação de diferenças finitas na equação que se deseja resolver numericamente. A equação algébrica (sistema linear) advinda deste processo de discretização é denominada de equação de diferenças finitas. Para que a resposta da equação de diferenças finitas (EDF) represente de forma satisfatória a solução da EDP original, é necessário que exigir que algumas propriedades sejam satisfeitas para que se possa de fato obter soluções numéricas limitadas (não infinitas) e uteis. A seguir, são discutidas estas propriedades de forma qualitativa e, na seção ^4.2.1, são apresentados critérios quantitativos para garantir respostas acuradas na formulação para a equação da onda adotada aqui.

Consistência

Uma EDF é consistente com uma EDP se a diferença entre a EDF e a EDP (o erro de truncamento da série de Taylor) vai a zero quando o espaçamento da malha e o incremento de tempo vão a zero independentemente, ou seja,

(36) $$|\mathrm{EDP}-\mathrm{EDF}|\to 0\mathit{}\mathrm{se}\mathit{}\mathrm{\Delta }t\to 0,h\to 0.$$

Pode acontecer da condição acima ser verdadeira somente quando uma certa relação entre o espaçamento da malha e o passo de tempo é satisfeita. Neste caso, o esquema é dito condicionalmente consistente. É claro que, se o erro de truncamento das aproximações da aproximação é conhecido e proporcional aos espaçamentos e incremento temporal (ou uma potência destes), como ocorre em todos os esquemas derivados com MDF, então a prova de consistência é direta, sendo os esquemas sempre consistente.

Estabilidade

Tal propriedade relaciona-se com o fato de que uma formulação numérica não pode amplificar artificialmente os erros numéricos da solução, o que degradaria completamente a solução numérica do problema ao longo do tempo. Neste sentido, se a solução de uma EDP é limitada, então a solução numérica deve ser limitada também. Como na maioria das vezes a solução de problemas físicos de interesse é limitada, pode-se definir estabilidade numérica como a seguir:

Um esquema de diferenças finitas é estável se produz uma solução limitada quando a solução exata é limitada, e, reversamente, é instável se produz uma solução ilimitada, sendo a solução exata limitada.

Há casos em que a estabilidade do esquema depende do espaçamento da malha e do incremento de tempo, sendo tais esquemas chamados de condicionalmente estáveis. Este é o caso do esquema numérico utilizado neste artigo e de grande parte dos esquemas de diferenças finitas. Na verdade, é o caso dos esquemas de diferenças finitas explícitos. Neste casos, é necessário que o espaçamento da malha e o incremento temporal estejam relacionados entre si de forma adequada, chamada de condição de estabilidade. O pior dos casos ocorre quando o esquema apresenta soluções numéricas ilimitadas para quaisquer valores de espaçamento da malha e incremento temporal, sendo em tais situações chamados de incondicionalmente instáveis. A chamada análise de estabilidade de von Neumann, uma aplicação da análise de Fourier, é o procedimento mais largamente empregado para provar a estabilidade de esquemas de diferenças finitas. Maiores detalhes sobre esta análise podem ser encontrados na Ref. [¹³].

Convergência

Isoladamente, é a propriedade mais importante. Ela exige que a solução numérica deve convergir para a solução do problema original (regido pela EDP) quando o espaçamento da malha e o incremento temporal tenderem a zero. Infelizmente, provar que um esquema é convergente é bem mais complicado do que provar a consistência e a estabilidade. Por este motivo, um teorema de importância prática imensa é conhecido como Teorema da Equivalência de Lax-Richtmyer [¹³, ¹⁴]. Tal teorema estabelece que um esquema é convergente se ele é estável e consistente. Assim, a tarefa de provar a convergência é substituída pelas tarefas mais simples de provar a consistência e a estabilidade.

4.2. Formulação explícita adotada

Nesta seção, será apresentada uma discretização simples da Eq. (10). Além disso, serão abordados critérios práticos para que tais soluções sejam adequadas.

O esquema abordado aqui utiliza expressões de quarta ordem no espaço (Eq. 35) e de segunda ordem no tempo (utilizando expressão análoga à Eq. 34 para o tempo). Em 2D, a equação de diferenças (que fornece a solução numérica) é obtida substituindo-se as aproximações e isolando no primeiro membro o termo do tempo mais avançando dado por $P_{i,k}^{n+1}$, obtendo-se:

(37) $$P_{i,j}^{n+1}={\delta }_{i,j}[-(P_{i+2,j}^n+P_{i,j+2}^n+P_{i-2,j}^n+P_{i,j-2}^n)$$ $$\mathrm{}+16(P_{i+1,j}^n+P_{i,j+1}^n+P_{i-1,j}^n+P_{i,j-1}^n)-60P_{i,j}^n]$$ $$\mathrm{}+2P_{i,j}^n-P_{i,j}^{n-1},$$

onde o parâmetro numérico definido acima é dado por

(38) $${\delta }_{i,j}=\frac{\mathrm{\Delta }{t}^2{C}_{i,j}^2}{12h}.$$

Um fato importante em relação à discretização no tempo é que esquemas em quarta ordem (utilizando expressão análoga à Eq. 35 para o tempo) são incondicionalmente instáveis [²⁰]. Por outro lado, as derivadas espaciais podem ser aproximadas em ordens mais elevadas sem problemas, o que é feito para minimizar a dispersão numérica (ver seção ^4.2.1).

Observa-se que, na formulação numérica apresentada, somente a velocidade c de propagação da onda é considerada, sendo este o único parâmetro que descreve o meio. Embora a densidade ρ apareça multiplicando a fonte, ela é considerada apenas no ponto onde ela é aplicada. Já a velocidade deve ser fornecida ponto a ponto no esquema, sendo este chamado de esquema heterogêneo (como definido na Ref. [²¹]).

A formulação de diferenças finitas definida pela Eq. 37 é uma formulação explícita, uma vez que a resposta pode ser calculada explicitamente nos tempos futuros em função de tempos anteriores. Dadas condições iniciais conhecidas (em n = 0), isto permite que a solução seja calculada de forma explicita no passo de tempo seguinte (em n = 1) em função de valores de $P_{i,k}^0$ conhecidos no tempo nulo e de valores $P_{i,k}^{-1}$ no passo anterior (todos iguais a zero, pois são tempos anteriores ao da aplicação da fonte, no qual o meio encontra-se em equilíbrio). De forma geral, conhecendo-se $P_{i,k}^n$ e $P_{i,k}^{n-1}$, pode-se calcular todos os valores de pressão no passo de tempo n + 1 e este processo continua iterativamente por N_t passos até que se atinja um determinado tempo final t_f = N_tΔt de análise. 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 se comparado à marcha explícita.

Os pontos de pressão discretos (sobre a malha de diferenças finitas) necessários para o cômputo de um ponto de pressão no tempo futuro definem o que se chama de estêncil do operador de diferenças finitas, conforme mostrado na Fig. 8.

Figura 8
Estêncil do esquema de diferenças finitas.

A seguir, será tratada a implementação computacional. Além da simplicidade mostrada acima para a obtenção de esquemas de diferenças finitas e da simplicidade conceitual do MDF, outra grande vantagem do método reside na simplicidade da programação de esquemas explícitos. Além disso, a solução numérica numérica é de boa qualidade, em especial para problemas de propagação de ondas, desde que sejam satisfeitas as condições discutidas na seção ^4.2.1, i.e., desde que se satisfaça o critério de estabilidade e se utilize espaçamentos da malha e incremento temporal pequeno o suficiente para que e se tenha uma boa precisão numérica.

No código computacional, em primeiro lugar, deve-se definir os parâmetros físicos associados ao problema. Estes são: o domínio de propagação da onda, ou seja, o tamanho X e Z do modelo 2D (em metros, por exemplo); as velocidades de propagação de onda no meio, definidas ponto a ponto no domínio; e a frequência de corte da fonte. Os parâmetros numéricos também precisam ser definidos de forma adequada, como será detalhado na seção ^4.2.1. Os principais parâmetros numéricos são: os espaçamentos da malha Δx e Δz, que define os números de pontos N_z×N_x, e o incremento temporal Δt, que define o número de passos N_t.

Um algoritmo básico é apresentado (Algoritmo 8) e explicado na sequência. Primeiro, deve-se declarar as variáveis necessárias para a solução do problema, definir os parâmetros físicos e numérico associados ou lê-los de arquivo. Caso sejam lidos de arquivo durante a execução do programa, as matrizes de N_z×N_x posições devem ser alocadas dinamicamente, após tais parâmetros terem sido lidos e alternativamente, de forma mais simples, tais matrizes podem ser declaradas se N_z e N_x são declarados como parâmetros. Na linha 8, o modelo de velocidade é lido para a matriz de velocidades. Na linha 8, deve-se realizar alguns cálculos iniciais, como: (1) zerar as matrizes P1 e P2 que irão guardar os valores de pressão na malha em instantes de tempo consecutivos, respectivamente em n−1 e n, necessários para a marcha explícita no tempo (cálculo de P3 em n + 1), (2) calcular parâmetros da fonte, como os dados pelas Eq. 12–14 e (3) a matriz δ é calculada (Eq. 38).

Algoritmo 1
Propagação de onda

A linha 8 dá início à marcha no tempo propriamente dita (laço sobre o número de passos de tempo n). Na linha 10, a fonte sísmica é aplicada no ponto kf×if da malha, utilizando-se uma função definida pelo usuário, com base na Eq. 11. A fonte é aplicada por meio da introdução do valor da fonte na posição correta da malha, por meio da matriz P1. Assim, ela será corretamente considerada quando executada a linha 14 na posição da fonte. Repare que a fonte só será aplicada se o passo de tempo for menor ou igual que a variável inteira N_fonte que especifica o número de passos de tempo onde a fonte é não nula, conforme definido pela Eq. 15. Deve-se utilizar N_fonte como o inteiro de valor mais próximo de t_M/Δt (t_M dado pela Eq. 16). Entre as linha 12–16, o operador de diferença finita, dado pela Eq. 37, é calculado em cada ponto da malha através de um laço duplo. Na linha 14, F(P2) representa todos os valores de P2 necessários para o cômputo de P3, dados pelos pontos entre colchetes na Eq.37 (Fig. 8).

Na linha 17, foram agrupadas em uma subrotina todo o tratamento das bordas, ou seja, a aplicação de segunda ordem nos pontos próximos ao contorno, a aplicação de p = 0 na superfície, das equações one-way da onda nas demais bordas, bem como as camadas de Cerjan. Detalhamento sobre esta subrotina é deixado como exercício para o leitor interessado em implementar o código utilizando uma linguagem de programação. As linhas 18–22 escrevem os snapshots em arquivos a cada N_snap passos, onde MOD(a,b) retorna o resto da divisão inteira a/b. As linhas 22–24 salvam o sismograma na matriz Sismo (cada um dos N_R receptores em (z,x) = i_prof,x_k), posteriormente escrita em arquivo. Na linha 25, os campos de pressão são atualizados para continuar a marcha. Dentro de um laço duplo sob toda a malha, utiliza-se as atribuições: P1(i, k)←P2(i, k), P2(i, k)←P3(i, k).

4.2.1. Critérios práticos para acurácia

Um esquema numérico pode satisfazer todas as propriedades discutidas acima (seção ^4.1.2) e ainda assim apresentar respostas numéricas com acurácia insuficiente para uma determinada aplicação. Intuitivamente, vê-se que a precisão da solução numérica depende dos parâmetros h e Δt adotados e, a princípio, quanto menores forem, melhor a resposta. Na realidade, o assunto é mais complexo como o leitor mais atento deve ter notado da discussão acima sobre estabilidade, em especial, do fato do nosso esquema ser condicionalmente estável.

Como dito, tais esquemas requerem a utilização de parâmetros h e Δt relacionados entre si, respeitando a condição de estabilidade. Entretanto, mesmo satisfazendo esta condição – que a rigor só garante que os erros numéricos não se tornem incontroláveis durante a marcha no tempo –, as soluções podem se desviar bastante da solução esperada, sendo necessário impor valores h e Δt ainda menores para garantir um determinado grau de acurácia da solução.

Um estudo detalhado relativo à obtenção do critério de estabilidade ou sobre a acurácia de formulações de diferenças finitas está fora do escopo deste trabalho. Ao leitor interessado em aprofundar-se no assunto, indica-se os mesmos livros de referência citados no início da seção ⁴. Aqui, apresenta-se, em lugar disso, alguns critérios práticos que, uma vez aplicados, garantem que as soluções numéricas vistas acima possuam estabilidade e precisão suficiente para inúmeras aplicações. Assim, frente aos parâmetros físicos do problema em questão (e.g., velocidade do meio, frequência da fonte), impõe-se a necessidade de os parâmetros numéricos satisfazerem os critérios apresentados abaixo.

Antes de apresentar os critérios, entretanto, faz-se necessário algumas explicações adicionais. As soluções numéricas de problemas ondulatórios estão sujeitas ao problema conhecido como dispersão numérica. Isto é, no caso de meios não dispersivos e isotrópicos, como os regidos pela equação acústica da onda, as velocidades de fase e de grupo (ver conceito na Ref. [⁸, ²²]) coincidem e são independentes da frequência e do ângulo de propagação. Quando discretiza-se o domínio do problema através do MDF, estas velocidades passam a ser função do espaçamento utilizado para a malha, além da velocidade de propagação do meio, da frequência e do ângulo de propagação. Para minimizar o nível de dispersão numérica a um patamar adequado (e imperceptível nas escala de tempo que será adotada) basta impor o seguinte critério:

(39) $$h\le \frac{c_{min}}{\alpha f_{corte}},$$

onde h é o espaçamento da malha, f_corte é a frequência de corte da fonte, isto é, a frequência máxima presente, c_min é a velocidade mínima presente no modelo de velocidades e α é um parâmetro que varia de acordo com a aproximação numérica utilizada. Assim definido, o parâmetro α representa o número de pontos discretos por comprimento de onda, onde se considera o menor comprimento de onda presente (referente à menor velocidade e à maior frequência).

Em esquemas de diferenças finitas utilizando operador de segunda ordem no tempo e quarta ordem no espaço, como o discutido aqui, valores de α = 5 são suficientes para um boa precisão. Utilizando operadores de segunda ordem no espaço, é necessário utilizar α = 10. Enfatiza-se, ainda, que estes valores dependem do esquema numérico considerado, podendo variar para outras formulações.

Para que a solução tenha estabilidade numérica, por outro lado, garantindo ainda um alto grau de precisão numérica, o seguinte critério deve ser satisfeito pelos parâmetros numéricos:

(40) $$\mathrm{\Delta }t\le \frac{h}{\beta c_{max}},$$

onde Δt é o intervalo de tempo, h é espaçamento da malha, c_max é a velocidade máxima presente no modelo e β é um parâmetro obtido de forma empírica.

A relação acima é similar ao critério de Courant-Friedrich-Levy (CFL). Ela expressa a necessidade do intervalo de tempo considerado ter de ser menor que o menor tempo que a onda leva pra atravessar uma célula de diferenças finitas (ou seja, considerando-se a maior velocidade no meio). Assim, β pode ser interpretado como o número de iterações necessárias para levar esta informação de um ponto a outro. Valores de β = 4 são suficientes para garantir a estabilidade da resposta numérica e um bom grau de precisão.

5. Aplicação ao “imageamento” sísmico

Esta seção discute como a equação acústica da onda e as soluções numéricas obtidas anteriormente podem ser utilizadas para construção de um método que é capaz de gerar imagens associadas a descontinuidades de propriedades físicas do meio, no caso em questão, à velocidade de propagação.

5.1. Aquisição sísmica

A seguir, apresenta-se uma aquisição sísmica sintética realizada com um programa escrito em Fortran 90 seguindo o algoritmo explicado anteriormente. O modelo rodado está representado na Fig. 9. Ele foi derivado do modelo desenvolvido pela Hess Corporation e disponibilizado por cortesia pela empresa, sendo mantido na Internet pela Society of Exploration Geophysicists (SEG) para que interessados possam utilizar. Trata-se de uma modelo offshore (em mar) interessante do ponto de vista geológico, uma vez que contém um corpo salino de alta velocidade (c = 4500 m/s) à esquerda e uma grande falha à direita. Originalmente ele é um modelo anisotrópico, como simetria vertical transversa (ver artigo II), embora seja utilizado aqui neste primeiro artigo da série considerando-se os parâmetros de anisotropia iguais a zero.

Figura 9
Modelo de velocidades acústico utilizado, baseado no modelo Hess.

O modelo Hess é composto originalmente por 3617×1501 pontos e uma lâmina d’água pouco profunda. O modelo utilizado aqui teve sua camada de água expandida e foi reamostrado para 2001×601 pontos, de modo a ser factível rodar os exemplos deste trabalho com a máquina disponível em tempo razoável. Como o espaçamento utilizado foi de h = 7 m, o tamanho total do modelo ficou em 14,0×4,2 km. A fonte empregada foi um pulso de Ricker de frequência f_corte = 40 Hz e o incremento temporal adotado foi de Δt = 3,8×10⁻⁴ s, sendo o número de passos no tempo adotado de N_t = 9000, totalizando um tempo de análise total de t_f = 3,42 s. Com isto, os parâmetros numéricos ficaram de acordo com o recomendado no tópico ^4.2.1, ou seja, o número α (pontos por comprimento de onda mínimo) não inferior a 5 e o número de passos de tempo β (iterações necessárias para levar a informação mais rápida de um ponto a outro da malha) não inferior a 4:

(41) $$\alpha ={\displaystyle \frac{c_{min}}{h{f}_{corte}}}={\displaystyle \frac{1500}{7\times 40}}\approx 5,4$$

(42) $$\beta ={\displaystyle \frac{h}{\mathrm{\Delta }t{c}_{max}}}={\displaystyle \frac{7}{3,8\times {10}^{-4}\times 4500}}\approx 4,1.$$

Tais valores garantem uma precisão numérica adequada. Como os valores estão próximos dos ideais estipulados, garante-se que a modelagem não utilize mais recursos computacionais que o necessário para o problema em questão. Por exemplo, rodar um modelo com as mesmas dimensões físicas e pelo mesmo tempo total de análise utilizando α = 10 (o que equivale a reduzir o h pela metade) implicaria, no caso 2D, em quadruplicar o número de pontos da malha, multiplicando o tempo de processamento pelo mesmo fator. Além disso, ao reduzir o h pela metade, o Δt também precisa ser reduzido pela metade para satisfazer a Eq. 40, implicando na necessidade de dar o dobro do número de passos de tempo. No total, o tempo de processamento teria aumentado 8 vezes, além do maior gasto de memória.

Na aquisição sintética rodada, um navio arrasta os receptores atrás de si, andando para a direita com a geometria de aquisição descrita a seguir. A cobertura do cabo de receptores é de 2,1 km, onde os receptores estão espaçados de 70 m (sendo 70 m o offset mínimo), totalizando 30 receptores localizados a esquerda da fonte sísmica (como mostrado esquematicamente na Fig. 9): a linha tracejada R representa os receptores e os asteriscos representam a fonte nas posições em que foram aplicadas. A distância entre os tiros foi também de 70 m e, no total, foram dados 165 tiros, com a posição inicial de x = 2,1 km. Somente as posições do primeiro tiro e dos tiros 10, 60, 110 e 160 são apresentadas na Fig. 9. Cada tiro demorou cerca de 7 minutos para roda em uma máquina com processador Intel i7 de 3.4 GHz.

Alguns instantâneos mostrando a distribuição do campo de pressão no domínio para um determinado tempo, para diferentes tempos, referente ao tiro 60, são apresentados na Fig. 10.

Figura 10
Snapshots da propagação da onda para o tiro 60 nos tempos (a) t = 1,33s, (b) t = 1,90s, (c) t = 2,47s e (d) t = 3,04s.

Na Figura 11, são apresentados os sismogramas dos tiros 10, 60, 110 e 160 (mostrados na Fig. 9), ou seja, as respostas ao longo do tempo (no eixo vertical) para cada receptor (no eixo horizontal), dispostos em ordem crescente de offset (distância entre fonte e receptor).

Figura 11
Sismogramas adquiridos nas posições indicadas na Fig. 9, referentes aos tiros de número (a) 10, (b) 60, (c) 110, (d) 160.

Através dos snapshots (instantâneos do campo de propagação em termos de pressão acústica), pode-se acompanhar a propagação da onda, desde a fonte até ser refletida e retornar aos receptores. Assim, pode-se identificar os diversos eventos coerentes presentes em um sismograma (que seja apresentado posteriormente). É fácil mostrar que as reflexões em interfaces planas são caracterizadas por hipérboles no sismograma, as chamadas hipérboles de reflexão. Quando os receptores não são planos, pode haver desvios da forma hiperbólica. Mesmo não sendo planos, eles podem ser aproximadamente planos localmente (como é o caso deste exemplo), dependendo do tamanho do cabo de receptores em comparação com a curvatura do refletor.

Em geral, sismogramas contém inúmeras hipérboles que podem ser facilmente identificadas em meio aos mais variados ruídos coerentes (que são eventos assim considerados uma vez que não se tratam de reflexões, como, por exemplo, as refrações). Neste exemplo, não foram inseridos ruídos adicionais (há, entretanto, ruídos numéricos, imperceptíveis nos instantâneos), e sismogramas representam as respostas puras das simulações de propagação, contendo todos os eventos associados a esta propagação, como por exemplo as múltiplas (eventos que ocorrem quando a onda atinge a superfície do modelo, sendo refletiva e retornando para baixo em direção aos refletores novamente, atingindo por fim os receptores mais uma vez).

O primeiro evento hiperbólico que aparece em cada um dos sismogramas sintéticos da Fig. 11 caracteriza a reflexão no fundo do mar, sendo as demais curvas associadas à reflexão nas demais interfaces e no sal, bem como reflexões múltiplas. A reta inclinada que aparece sempre na mesma posição em todos os sismogramas é a onda direta que se propaga na água, desde a fonte até os receptores. Quanto mais longe o receptor está da fonte (maior o offset), mais tempo a onda irá levar para chegar ao receptor. No tiro 60 (Fig. 11(b)), como exemplo de identificação de eventos nos sismograma sintético, chama-se a atenção para o sexto evento de reflexão, com amplitude acentuada, aproximadamente entre os tempos 2,25 e 2,5 s. Pode-se ver claramente através dos snapshots (em especial, da Fig. 10(a)) que ele caracteriza a reflexão no corpo de sal.

5.2. Mapeamentos dos refletores

Nesta seção, dentro do escopo do método sísmico, serão apresentadas as ideais fundamentais sobre o “imageamento” (mapeamento dos refletores) utilizando a equação da onda aplicada no sentido inverso do tempo, denominado de Migração Reversa no Tempo pré-empilhamento na literatura especializada da área de Geofísica Aplicada. Posteriormente, são descritas em detalhes as etapas do algoritmo de mapeamento e, na sequência, as alterações necessárias no algoritmo de modelagem direta para a sua programação em computador utilizando uma linguagem de programação de alto nível (neste trabalho, foi utilizada a linguagem Fortran, como já mencionado).

Neste artigo, diferente do procedimento padrão realizado rotineiramente na indústria, o algoritmo será aplicado a dados sintéticos, gerados com o algoritmo de modelagem direta. Enfatiza-se que este tipo de migração (sintética) é muito utilizado para verifica a eficácia de um determinado método, por exemplo, frente a uma configuração geológica de interesse.

Matematicamente, o problema do imageamento é um problema inverso, em contraposição ao problema direto da propagação da onda. No problema direto, as propriedade do meio são conhecidas e se deseja obter as respostas relativas à propagação da onda neste meio. No problema inverso, a resposta da propagação da onda no meio é conhecida e queremos utilizá-la, de forma inversa, para obter informações sobre o meio.

Contudo, o problema do imageamento dentro do escopo do método sísmico abordado neste trabalho é um problema inverso aproximado, uma vez que ele utiliza o que se chama de macromodelo de velocidades como informação a priori (input do algoritmo). Como discutido anteriormente na seção ², o objetivo principal dos diferentes métodos aplicados para o imageamento em sísmica de reflexão, conhecidos genericamente como métodos de migração, é o reposicionamento da energia sísmica medida para a suas corretas posições em subsuperfície, retirando as distorções introduzidas pela propagação da onda. Ou seja, objetiva-se obter uma imagem dos refletores e difratores em subsuperfície, relacionados com a geologia local.

Dentro do longo fluxo de processamento utilizado rotineiramente pela indústria nestes problemas, é realizado previamente ao imageamento um procedimento chamado de análise de velocidades, onde se objetiva estimar o macromodelo que será utilizado na etapa de migração. Ele poderá ser utilizado para a migração ou, o que é mais comum, servirá de base para algum procedimento mais avançado para geração do macromodelo final, e.g., tomografia ou inversão da forma de onda completa (full waveform inversion, um procedimento de inversão não linear). É importante mencionar que a imagem gerada com a utilização do macromodelo, ou seja, o resultado da migração, fornece indícios sobre a qualidade do macromodelo. Quando o “foco” das imagens é ruim ou existem refletores se cruzando, por exemplo, sabe-se que há problemas. Quando isto ocorre, é comum que o processo de imageamento, quando aplicado em campo, seja feito mais de uma vez, até se chegar a uma solução considerada de qualidade suficiente pelo interprete (o geofísico responsável pela interpretação da imagem com o objetivo final de decidir a locação de um poço exploratório que irá comprovar ou não a existência de hidrocarbonetos no alvo previamente selecionado).

Existem muitos métodos de migração diferentes, sendo um dos mais populares a chamada migração Kirchhoff, uma técnica que utiliza uma solução integral do problema de propagação de ondas para imagear estruturas. Ela é muito utilizada porque é uma técnica extremamente rápida onde as soluções integrais são implementadas de forma muito eficiente através de somatórios, sem marcha no tempo. Esta técnica, entretanto, não é adequada para utilização em modelos com algo grau de complexidade, como aqueles associados a regiões envolvendo tectônica salina, por exemplo. Nestes casos, em que o modelo de velocidade é complexo, os métodos de migração sísmica baseados na equação da onda vêm ganhando destaque. Não somente os crescentes desafios exploratórios, mas também o desenvolvimento de computadores mais rápidos (clusters) e com maior capacidade de armazenamento, tonaram este um ramo importante. Estes métodos são conhecidos genericamente por métodos de migração baseados na equação da onda. Nestes métodos, os algoritmos de modelagem e de migração são utilizados em conjunto, apresentando as mesmas características.

O que eles têm em comum é o fato de resolverem alguma aproximação da equação diferencial da onda (ou mesmo a equação completa), através da propagação direta do campo, diferentemente do método de migração Kirchhoff, que é baseado em soluções integrais desta equação, nas quais o campo de onda não é propagado. Tais métodos podem diferir quanto ao domínio utilizado (tempo ou frequência), quanto ao tipo de equação da onda utilizada (acústica ou elástica) e quanto à aproximação utilizada para resolver numericamente a equação da onda, como por exemplo a aproximação unidirecional que considera apenas o campo de onda descendente.

Neste trabalho, será utilizada a equação acústica completa da onda, resolvida no domínio do tempo, o que caracteriza a Migração Reversa no Tempo. Especificamente, se utilizará a técnica que migra tiros separados e empilha as imagens resultantes de cada tiro ao final para obter a imagem completa, ou seja, a técnica pré-empilhamento. A migração pré-empilhamento apresenta resultados muito melhores do que a técnica pós-empilhamento, onde os dados são empilhados antes da realização da migração. O processo de empilhamento antes da migração sofre de problemas sérios para “focar” a energia em regiões complexas, o que terminam por degradar a imagem final. Entretanto, como a técnica de migração pós-empilhamento é aplicada apenas uma vez (a toda a seção empilhada), ele demanda muitos menos tempo de computação que a migração pré-empilhamento, sendo escolhida em alguns casos práticos por esta vantagem.

A seguir, serão abordadas as etapas do método de migração utilizando propagação do campo de ondas reversamente no tempo, ou seja, a Migração Reversa no Tempo pré-empilhamento que é a técnica que será aplicada neste artigo.

Pode-se dizer que o método de imageamento discutido neste artigo é o mais simples e intuitivo do ponto de vista conceitual para o imageamento sísmico.

Por simplicidade, aqui o método será discutido levando em conta a aplicação a um único tiro (sismograma de tiro), com ênfase nas ideias fundamentais. Para obter uma imagem completa da região de interesse, é necessário que se aplique o algoritmo a diversos tiros e posteriormente as imagens são somadas (processo chamado de empilhamento).

Em primeiro lugar, como dito anteriormente, o algoritmo desenvolvido aqui será aplicado a um exemplo sintético, ou seja, os sismogramas que serão utilizados para imagear estruturas em subsuperfície vêm de uma modelagem numérica do problema direto. Isto evita a necessidade de pré-processar os dados antes de migrar, de forma a ser possível se concentrar na parte conceitual do método e ainda assim chegar a um resultado final de excelente qualidade.

Tudo o que será explicado aqui se aplicada da mesma forma a dados reais. De uma forma ou de outra, a migração pressupõe a existência de um dado a ser migrado (no nosso caso, os sismogramas de tiro, como os mostrados no tópico ^5.1). Além disso, é necessário o conhecimento prévio do macromodelo de velocidades, como discutido acima. Como irá se trabalhar com dados sintéticos, o modelo de velocidade é conhecido. Ao invés do modelo de velocidades conhecido (que possui descontinuidades de propriedades física), serão utilizados modelos suavizados, onde as reflexões serão bastante atenuadas. Isto é feito para evitar distorções nos tempos de trânsito, devendo ser suavizado o modelo de vagarosidade (inverso da velocidade), como discutido na Ref. [¹²].

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 que se deseja imagear (seja um sismograma real, vindo de campo, ou um sintético) é 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 final de análise até o tempo inicial). Uma imagem para esta segunda etapa é a seguinte: considerando-se a onda que foi propagada da fonte até os receptores com um filme, a segunda etapa nada mais é do que voltar o filme para trás, até o início, de forma que os instantâneos vão ficando com o formato da fonte sísmica no exato local em que foi aplicada. Aplica-se, então, uma condição de imagem responsável por fazer com que a energia propagada no sentido inverso do tempo pare no local em que foi gerada através de uma reflexão ou difração. Em outras palavras, a condição de imagem é responsável por posicionar a energia na sua correta posição em subsuperfície.

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). Esta condição está representada esquematicamente na Fig. 12. Ela é chamada de condição e imagem por tempo de excitação.

Figura 12
Condição de imagem.

A seguir, detalha-se cada uma das etapas citadas acima, com ênfase na sua aplicação computacional, sendo os respectivos algoritmos discutidos em detalhes. Posteriormente, com o objetivo de exemplificar e tornar claro todo o processo, os algoritmos serão aplicados aos dados sintéticos apresentados no tópico ^5.1.

5.2.1. Implementação computacional

Cada uma das etapas mencionadas acima é implementada computacionalmente utilizando o código de diferenças finitas discutido na seção ^4.2 com algumas modificações pontuais adequadas para a obtenção de matriz de tempo de trânsito (MTT) e da imagem, respectivamente, na primeira e na segunda etapa, como será detalhado a seguir.

1a Etapa: propagação no sentido direto

Nesta etapa, objetiva-se obter os tempos que a onda direta demora para chegar a cada ponto do domínio, ou seja, a chamada matriz de tempo de trânsito (MTT). Para tal, o algoritmo de propagação de ondas é utilizado. Porém, é necessária lançar mão de algum critério para identificação da onda direta no algoritmo.

Um critério simples que pode ser utilizado com bons resultados, embora tenha os seus problemas inerentes, é o critério da amplitude máxima. Com este critério, considera-se o tempo de trânsito da onda com maior amplitude que passa em um determinado ponto do modelo durante o tempo de aquisição. Ou seja, a onda com maior amplitude é interpretada como a onda direta, embora se saiba que há situações em que isto não é verdade. Entretanto, este critério funciona bastante bem porque a amplitude da onda decai a medida que a onda se propaga (devido a divergência esférica) e, além disso, ao sofrer reflexões e refrações, a onda perde energia e, em geral, tem a sua amplitude diminuída. Com isto, a onda direta costuma ser a onda de maior amplitude se comparada às ondas que irão passar posteriormente no mesmo ponto de observação.

No Algoritmo 2, é apresentado o pseudo-código para obtenção da MTT. A primeira modificação que deve ser feita no algoritmo original (Algoritmo 1), conforme mostrado na linha 5, se deve à utilização do modelo de velocidades suavizado, para evitar a interferência de ondas refletidas em determinados pontos do domínio que interfiram construtivamente e gerem uma onda com amplitude grande que pode ser erroneamente avaliada como uma onda direta.

Figura 13
Matriz de tempo de trânsito utilizando o modelo de velocidades (a) não suavizado e (b) suavizado.

Na Figura 13 são mostradas, comparativamente, matrizes de tempo de trânsito com condição de máxima amplitude, com e sem a utilização do modelo de velocidades suavizado. Como pode ser visto, a MTT que utilizou o modelo de velocidades suavizado apresenta uma continuidade maior, caracterizando sua maior qualidade em comparação a que usa o modelo não suavizado.

Algoritmo 2
1a Etapa do imageamento

Na migração, não se deseja obter sismogramas (e, em geral, nem instantâneos da propagação). Assim, as linhas correspondentes a estas tarefas são substituídas por um código para obter a MTT, mostrado entre as linhas 17–24. Repare que um matriz Amp(N_z,N_x) auxiliar é utilizada para armazenar as amplitudes atuais dos tempos registrados para cada posição da malha. Tais amplitudes são verificadas em cada ponto da malha (a cada passo de tempo) e, caso a amplitude atual no ponto seja maior que a armazenada, registra-se então o novo tempo associado a esta amplitude. Quando a maior amplitude passa pelo ponto, o tempo correspondente é armazenado. O processo deixa então de registrar novos tempos, uma vez que eles estarão associados a amplitudes menores. Ao final da marcha no tempo, a MTT é escrita em arquivo de forma a poder ser utilizada na etapa seguinte da migração.

Mais uma vez, é importante lembrar que o procedimento descrito deve ser aplicado a cada um dos tiros da migração pré-empilhamento.

2a Etapa: propagação no sentido inverso

A seguir, a segunda parte da migração é apresentada no Algoritmo 3. Como na etapa anterior, um modelo de velocidades suavizado deve ser utilizado. Agora o objetivo é fazer com que as reflexões e difrações associadas à energia propagada reversamente no tempo sejam atenuadas, evitando que a condição de imagem produza artefatos advindos destas ondas.

Nesta etapa, adicionalmente, é necessário ler os arquivos contendo a matriz de tempo de trânsito e sismograma (armazenadas nas matrizes MTT e SIS), relativo ao tiro a ser migrado, como mostrado na linha 7. Além disso, o laço da marcha no tempo é realizado no sentido de tempo inverso, do tempo final para o inicial, conforme apresentado na linha 8. O sismograma deve ser aplicada como fonte, como discutido anteriormente, para que a energia registrada seja reintroduzida no domínio. Assim, a cada passo de tempo, aplica-se o sismograma (armazenado na matriz SIS(N_t,N_x)), na posição em que foram medidos, referentes a cada um dos traços sísmicos, ou seja, na profundidade i_prof e posição x₁ até x_k, com k = 1 até k = N_R, como feito no laço entre as linhas 9–11.

Algoritmo 3
2a Etapa do imageamento>

O laço duplo, varrendo toda malha, entre as linhas 18–23 é a condição de imagem. Como se vê, a cada passo de tempo, cada posição da malha é testada para verificar se o tempo (no caso dado pelo índice n de passo no tempo) é igual a MTT naquele ponto. Em caso afirmativo, a amplitude propagada inversamente no tempo (ou seja, P3) é salva na matriz imagem, na posição correspondente. Ao final da marcha no tempo, a matriz imagem é escrita em arquivo. Repare que esta imagem é referente a um único tiro, sendo, portanto, uma imagem parcial, capaz que imagear os refletores próximos à região do tiro. Da mesma forma, todos os demais tiros precisam ser migrados, devendo todas elas serem somadas para obter a imagem final.

5.3. Resultados do imageamento

Nesta seção, apresenta-se um exemplo completo de imageamento utilizando algoritmos escritos em Fortran 90 seguindo a lógica discutida acima. Como mencionado anteriormente, o método de imageamento utilizado é conhecido como Migração Reversa no Tempo pré-empilhamento. Neste exemplo serão migrados os 165 tiros do exemplo de aquisição apresentado no tópico ^5.1. Por consistência, a mesma geometria de aquisição e os mesmos parâmetros são utilizados.

Na Figura 14 são mostradas as imagens dos tiros 10, 60, 110 e 160, colocando-se o modelo de velocidades como fundo (para comparação). Como dito anteriormente, cada imagem é capaz de recuperar apenas as informações próximas aos receptores, sendo imagens parciais de toda a área em estudo. Somente após o empilhamento de todos o tiros é que obtém-se uma imagem satisfatória de toda seção. O processo de empilhamento é responsável também por melhorar a razão sinal-ruído da imagem final em relação às imagens individuais de cada trecho. Isto se deve ao fato dos tiros serem dados próximos uns aos outros, de forma a haver eventos redundantes.

Figura 14
Imagens de tiro único para cada um dos tiros mostrados na Fig. 9, ou seja, os tiros (a) 10, (b) 60, (c) 110 e (d) 160.

Assim, os eventos associados a reflexões e difrações de fato são amplificados, uma vez que tendem a aparecer em todas as imagens próximas, enquanto os artefatos tendem a se cancelar. O resultado final da migração dos 165, após o empilhamento de todos os tiros, é mostrado na Fig 15. Enfatizamos que a imagem está apresentada como saída do algoritmo de empilhamento, no qual as matrizes de imagem dos diferentes tiros foram simplesmente somada e nenhum filtro para remover ruídos foi aplicado. A migração dos 165 tiros levou cerca de 20 horas no total, sendo que a primeira e a segunda etapas foram rodadas em paralelo em um processador multinúcleo de 3.4 GHz.

Figura 15
Imagem final do modelo Hess após o empilhamento dos 165 tiros.

6. Conclusões

Neste trabalho, foi apresentada de forma didática e detalhada a parte teórica e a construção de um algoritmo de propagação de ondas acústicas utilizando o Método das Diferenças Finitas. Foi discutida também a sua aplicação ao mapeamento de refletores utilizando propagação de ondas reversas no tempo, sendo tal método conhecido como Migração Reversa no Tempo pré-empilhamento na literatura especializada de Geofísica Aplicada. Trata-se de um dos métodos de imageamento mais poderosos aplicados pela indústria do petróleo, que vem sendo utilizado para imagear regiões extremamente desafiadoras como, e.g., o pré-sal.

Todas as ideias fundamentais foram discutidas em detalhes, tendo os algoritmos sido apresentados, permitindo ao leitor uma visão detalhada de um dos principais métodos utilizados na exploração de hidrocarbonetos atualmente. Diversos exemplos e figuras visaram ilustrar ideias de compreensão mais difícil, mesmo para o público da área de Geofísica. A partir do exposto, espera-se que qualquer pessoa com conhecimentos de programação e física seja capaz de compreender e implementar os algoritmos apresentados e reproduzir os exemplos, preenchendo-se uma lacuna de material sobre o tema para um público de físicos, em especial aqueles que almejam especializar-se em Geofísica, um quantitativo cada vez maior de jovens, devido às perspectivas de emprego em empresas.

Além disso, a simplicidade da resolução da equação acústica da onda pelo MDF permite a abordagem do tema em disciplinas da graduação. Sua aplicação ao imageamento constitui-se em um desafio interessante para o estudante, uma vez que articula conhecimentos computacionais, de física pura e aplicada em um único problema, enquanto introduz o estudante a um tema de ponta em uma área correlata da Física.

O imageamento em meios mais complexos, que demandam a utilização de algoritmos de propagação mais elaborados – tanto no que se refere a equações mais refinadas, e.g., a equação elástica da onda com anisotropia, como outros algoritmos utilizando MDF ou outros métodos numéricos –, pode ser abordado utilizando-se as mesmas ideias básicas discutidas neste trabalho. Assim, espera-se que ele possa servir de motivação para o aprendizado e aplicação em casos mais complexos. Pode ainda servir de motivação para aplicações em casos distintos como, por exemplo, o imageamento com ultrassom em medicina ou a realização de ensaios não destrutivos em engenharia.

Agradecimentos

Agradeço aos Drs. Cleberson Dors e Cid S. G. Monteiro, autores do código de visualização utilizados e aos Drs. Jaime F. Vilas da Rocha e Luis Juracy R. Lemos pelas leituras criticas do manuscrito.

References

• ^1.
  E. Baysal, D.D. Kosloff e J.W.C. Sherwood, Geophysics 48 , 1514 (1983).
• ^2.
  G.A. Mcmechan, Geophysical Prospecting 31 , 413 (1983).
• ^3.
  N.D. Whitmore, Iterative depth migration by backward time propagation , disponível em: https://doi.org/10.1190/1.1893867
  » https://doi.org/10.1190/1.1893867
• ^4.
  P. Kearey, M.l Brooks e I. Hill, Geofísica de exploração (Oficina de Textos, São Paulo, 2009).
• ^5.
  O. Yilmaz, Seismic Data Analysis. Processing, Inversion, and Interpretation of Seismica Data (Society of Exploration Geophysicists, Tulsa, 2008).
• ^6.
  A.V. Andrade-Neto e D.D. Andrade, Revista Brasileira de Ensino de Física 43 , e20200350 (2021).
• ^7.
  C.P.A. Wapenaar e A.J. Berkhout, Elastic Wave Field Extrapolation: Redatuming of Single- and Multi-Component Seismic Data (Elsevier Science, New York, 1989), v. 2.
• ^8.
  M. Nussenzweig, Curso de Física Básica (Edgard Blücher, São Paulo, 1996), v. 2, 3ª ed.
• ^9.
  L. Di Bartolo, C. Dors e W.J. Mansur, Geophysics 77 , T187 (2012).
• ^10.
  C. Cerjan, D. Kosloff, R. Kosloff e M. Reshef, Geophysics 50 , 705 (1985).
• ^11.
  J.P. Berenger, Journal of computational physics 114 , 185 (1994).
• ^12.
  D. Loewenthal, P.L. Stoffa e E.L. Faria, Geophysics 52 , 1007 (1987).
• ^13.
  J.W. Thomas, Numerical partial differential equations: finite difference methods (Springer Science & Business Media, Berlim, 1995), v. 22.
• ^14.
  J.C. Strikwerda, Finite difference schemes and partial differential equations (SIAM, Pensilvânia, 2004), 2ª ed.
• ^15.
  R.J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems (SIAM, Pensilvânia, 2007).
• ^16.
  I.P. Mysovskih, Lectures on Numerical Methods (Wolters-Noordhoff Plublishin Groningen, Netherlands, 1969).
• ^17.
  V.S. Ryaben’Kii e S.V. Tsynkov, A theoretical introduction to numerical analysis (Chapman & Hall/CRC, New York, 2007).
• ^18.
  E. Isaacson e H.B. Keller, Analysis of numerical methods (Dover Inc, Illinois, 1994), 1ª ed.
• ^19.
  K.J. Bathe, Finite Element Procedures (Prentice-Hall, New Jersey, 1996), 2ª ed.
• ^20.
  L. Anne, P. Joly e Q.H. Tran, Computacional Geociences 4 , 207 (2000).
• ^21.
  K.R. Kelly, R.W. Ward, S. Treitel e R.M. Alford, Geophysics 41 , 2 (1976).
• ^22.
  K.F. Graff, Wave Motion in Elastic Solids (Dover Publication, Mineola, 1991).

• 1
  n-ésima ordem na variável x significa que o erro de truncamento da série de Taylor da aproximação via DF seja da ordem da n-ésima potência do equivalente espaçamento da malha, isto é, Δxⁿ.

Datas de Publicação

• Publicação nesta coleção
  06 Ago 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  02 Jun 2021
• Revisado
  13 Jul 2021
• Aceito
  13 Jul 2021