Modelagens direta e inversa de dados mCSEM 3D

Modelagens direta e inversa de dados mCSEM 3D

Frayzer Lima de Almeida^I; Luiz Rijo^II

^IUniversidade Federal do Pará, Campus Universitário de Castanhal, Faculdade de Matemática, Av. Universitária, s/n, Jaderlândia, 68745-000 Castanhal, PA, Brasil. Tel./Fax: [phone] - E-mail: [e-mail]

^II in memoriam

RESUMO

Marine Controlled Source Electromagnetic - mCSEM é um método geofísico eletromagnético que nos últimos dez anos vem sendo usado na prospecção de hidrocarbonetos com bastante êxito. Este método consiste em um dipolo elétrico horizontal (DEH) localizado um pouco acima do assoalho marinho, operando em baixa frequência (0,1-1,0 Hz) e receptores regularmente distribuídos no fundo do mar que captam os campos eletromagnéticos provenientes da difusão de energia gerada pelo dipolo transmissor. Neste trabalho vamos apresentar o problema direto do método mCSEM 3D, propondo soluções numéricas, através do método dos elementos finitos tridimensionais, para modelos geoelétricos mCSEM 3D. Para fins de análise de coerência, os resultados obtidos são comparados com soluções disponíveis na literatura. Em seguida, apresentaremos a inversão de um de seus modelos segundo uma proposta de metodologia de inversão juntamente com a proposta de solução direta para o mCSEM 3D, acima mencionada, realizando assim a inversão de um modelo geoelétrico do mCSEM 3D para duas frequências.

Palavras-chave: geofísica aplicada, métodos eletromagnéticos, elementos finitos, inversão geofísica, método de Marquardt, Método Eletromagnético Marinho de Fonte Controlada (mCSEM), exploração de petróleo, águas profundas.

ABSTRACT

Marine Controlled Source Electromagnetic - mCSEM is a new geophysical technology which has been successfully applied to hydrocarbon exploration in the last ten years. This method consists of a horizontal electric dipole (HED) moving just above the sea floor, operating at low frequencies (0.1-1.0 Hz). The fixed receivers distributed regularly on the sea floor capture the electromagnetic fields related with the diffusion of energy generated by the dipole transmitter, the sea water and the sediments, including the petroleum reservoir. In this work we present a finite element algorithm for 3D mCSEM modeling. The performance of the algorithmis illustrated by means of some selected examples in literature. In this work we propose also an algorithm for inversion of three-dimensional mCSEM synthetic datafor two frequencies. Results that illustrate the effectiveness of the proposed algorithm are also discussed.

Keywords: applied geophysics, electromagnetic methods, finite elements, geophysical inversion, Marquardt method, Marine Controlled Source Electromagnetic (mCSEM), oil exploration, deep waters.

INTRODUÇÃO

A prospecção de hidrocarbonetos em águas profundas apresenta risco exploratório muito elevado. Em decorrência disto, faz-se necessário o desenvolvimento e aprimoramento de novos métodos não-sísmicos para auxiliar a sísmica, este sendo reconhecidamente o principal método geofísico na exploração de petróleo. Neste sentido, o método mCSEM é bastante promissor (Constable & Srnka, 2007).

Na inversão de dados mCSEM com modelos 1D Christensen & Dodds (2007), que solucionam o problema inverso usando uma variação do método numérico não linear Gauss-Newton (Tarantola & Valette, 1982) juntamente com estabilizadores. Neste trabalho os dados observados sintéticos usados na inversão são as amplitudes e as fases das componentes eletromagnéticas E_x e H_y, nas frequências de 0,05 Hz, 0,15 Hz e 0,25 Hz, com contaminação por ruídos, sendo 2% para as amplitudes dos campos E_x e H_y e 3,5% para suasfases. Constable & Weiss (2006) que solucionam o problema inverso usando o método de OCCAM. Neste caso, os dados observados usados na inversão são as amplitudes da componente E_f com 10% de ruído, as amplitudes de E_r também com 10% de ruído e finalmente as resistividades aparentes e as fases (MT) com 10% de ruído e 2,9% de ruído, respectivamente. Para conhecer desenvolvimentos sobre o método de OCCAM sugerimos o trabalho de Constable et al. (1987).

Para a reconstituição de modelos do mCSEM 2D tem-se o trabalho de MacGregor et al. (2001), que usa para inversão oalgoritmo de OCCAM.

Para modelagem direta dos dados de mCSEM 3D tem-se os trabalhos de Ueda & Zhdanov (2005), que solucionam o problema direto através do método de equação integral, e o trabalho de Souza (2007), que soluciona o problema direto através do método dos elementos finitos. Este último soluciona o sistema linear associado através do método do gradiente biconjugado (Souza et al., 2005). Para a inversão de dados mCSEM 3D tem-se o trabalho de Zhdanov & Yoshioka (2005), que possui o problema direto solucionado através de equação integral e usa, para solucionar o problema inverso, uma variação do método não linear do gradiente conjugado juntamente com pesos regularizadores, ambos calculados através de computação paralela. Os dados observados usados na inversão são as componentesE_x, E_y e E_z do campo elétrico nas frequências de 0,25 Hz, 0,5 Hz, 1,0 Hz, 2,0 Hz e 4,0 Hz.

Primeiramente, para a modelagem direta do mCSEM 3D, este trabalho contribui com uma curta modificação na formulação proposta por Souza (2007), através do método dos elementos finitos tridimensionais usando os potenciais magnético A e elétrico Φ. Para fim de avaliação numérica dos nossos resultados de modelagem direta do mCSEM 3D proposta, apresentaremos comparações entre as amplitudes normalizadas da componente E_x (in-line ) (Souza, 2007), provenientes de modelos geoelétricos do mCSEM 1D, 2D e 3D. Secundariamente, este trabalho contribui com a proposta de um algoritmo parainversão de dados sintéticos do mCSEM 3D segundo o método de Marquardt (Press et al., 1992), usando apenas a componente E_x do campo elétrico, para uma frequência, localizadono assoalho marinho, juntamente com a modelagem direta do mCSEM 3D também proposta inicialmente. Diferentemente daliteratura acima citada, em que esta usa para a modelagem direta do mCSEM 3D o método de equação integral e uma variação do método não linear do gradiente conjugado para a inversão de dados do mCSEM 3D através das componentes E_x, E_y e E_zdo campo elétrico, para cinco frequências.

Este trabalho está organizado da seguinte forma: inicialmente, na seção de características preliminares do mCSEM 2D e 3D, são expostas noções, através de figuras e descrições, sobre os ambientes geoelétricos do mCSEM bidimensional e tridimensional, assim como a citação dos métodos numéricos aplicados para suas modelagens diretas. Em seguida, na seção de metodologias, apresentamos o problema direto do mCSEM 3D solucionado pelo método dos elementos finitos, usando os potenciais EM, e o problema inverso do mCSEM 3D solucionado pelo método de Marquardt (Press et al., 1992), juntamente com a introdução de vínculos absolutos. A seguir, na seção de resultados, são expostas e comparadas as soluções diretas do mCSEM 1D, 2D e 3D, sendo esta última proposta por este trabalho, assim como a apresentação da reconstituição de um modelo geoelétrico tridimensional, para duas frequências, usando o algoritmo de inversão de dados sintéticos do mCSEM 3D também proposto. Em seguida, na seção de conclusões, tanto as contribuições, quanto as análises dos resultados obtidos são evidenciadas deforma resumida.

CARACTERÍSTICAS PRELIMINARES DO mCSEM 2D e 3D

Esta seção tem a finalidade de apresentar noções preliminaresdos ambientes geoelétricos referentes às características geológicas brasileiras de prospecção de hidrocarbonetos em águas profundas do mCSEM 2D e 3D. Também se destina a apresentar uma breve descrição metodológica acerca da modelagem do mCSEM 2D e 3D no sentido de como solucionar os campos eletromagnéticos do referido método.

Características da modelagem do mCSEM 2D

O modelo mCSEM bidimensional (Fig. 1) é formado por camadas estratificadas. A camada do ar possui resistividade ρ₀ infinita. A camada do mar, uma resistividade ρ₁ em torno de 0,3 Ωm e espessura h₁ que varia entre 1 a 3,5 km. A camada das rochas sedimentares encaixantes possui resistividade ρ₂ em torno de 1,0 Ωm. Uma camada plano paralela representa um reservatório de hidrocarbonetos que possui resistividade ρ₃ entre 10 e 100 Ωm, com uma espessura entre 50 e 100 m e limitada na direção do eixo x. A fonte é um dipolo elétrico horizontal(DEH) em uma profundidade h₀ com seu eixo na direção x e, finalmente, os receptores são representados por asteriscos.

Os campos eletromagnéticos, mCSEM 2D, referentes ao modelo acima descrito, são solucionados numericamente a partirdas equações de Maxwell no domínio da frequência pelo método dos elementos finitos bidimensionais (Rijo, 2004).

Características da modelagem do mCSEM 3D

O modelo tridimensional (Fig. 2) é novamente formado por camadas estratificadas. A camada do ar possui resistividade elétrica infinita. A camada do mar uma resistividade em torno de0,3 Ωm. A camada das rochas sedimentares encaixantes possui resistividade em torno de 1,0 Ωm. Um volume tabular, paralelo às camadas estratificadas, representando um reservatório de hidrocarbonetos, possui resistividade entre 10 e 100 Ωm comuma espessura entre 50 e 100 m e limitado nas direções dos eixos x, y e z. A fonte também é um dipolo elétrico horizontal (DEH) com seu eixo na direção x, a uma altura fixa do assoalho marinho, podendo ser entre 20 e 50 m e podendo também ter várias localizações em relação aos eixos x e y. Finalmente o modelo é completado pelos os receptores localizados e afixados no assoalho marinho, que registram os campos eletromagnéticos. Geralmente em cada linha de receptores, que são paralelas ao eixo x, os mesmo podem apresentar espaçamentos de 500 m.

Estes campos eletromagnéticos, mCSEM 3D, são solucionados numericamente a partir das equações de Maxwell no domínio da frequência pelo método dos elementos finitos tridimensionais (Souza, 2007), como também pelo método de equação integral tridimensional (Ueda & Zhdanov, 2005).

METODOLOGIAS

O problema direto do mCSEM 3D segundo os potenciais EM (A,Φ)

Usando as equações de Maxwell no SI em um meio não magnético, no domínio da frequência, em termos dos campos secundários elétrico E^s e magnético H^s, temos:

em que ε₀ = 8,854 · 10^-12 F/m e µ₀ = 4π· 10^-7 H/m são a permissividade elétrica no vácuo e a permeabilidade magnética no vácuo, respectivamente. σ^p e Δσ são a condutividade elétrica do meio homogêneo e a diferença de condutividade entre o meio homogêneo e a heterogeneidade, respectivamente. J^s = ΔσE^s é a densidade de correntes secundárias e w = 2πf a frequência angular, enquanto f a frequência linear. Em virtude das Equações (3) e (4) podemos expressar os campos E^s e H^s em termos do potenciais secundários magnético A^s e elétrico Φ^s = iωµ₀Φ^s da seguinte forma:

Das Equações (5) e (6), a Equação (2) torna-se a equaçãorotacional do rotacional

A Equação (7), quando expressa segundo o método dos elementos finitos, resulta em uma matriz assimétrica que pode apresentar modos espúrios e instáveis numericamente (Souza, 2007). Para evitar esta dificuldade, aplica-se a identidade vetorial

juntamente com a incorporação do calibre de Coulomb:

dando origem à equação

para a qual a sua forma discretizada é estável numericamente (Souza, 2007). Na derivação da Equação (7) definiu-se o potencial escalar reduzido F^s = iwm₀f^s, o que resulta em uma matriz de elementos finitos simétrica. Tomando o divergente da Equação (7) e considerando a propriedade ∇·∇ × H = 0, tem-se

Desta forma as Equações (8) e (9) são solucionadas simultaneamente através do método dos elementos finitos (Souza,2007), do qual se obtém um sistema linear associado a uma matriz complexa, simétrica e esparsa. Este sistema é resolvido numericamente pelo método do gradiente biconjugado (Souzaet al., 2005), resultando assim nos potenciais A^s e Φ^s.

Na formulação do problema direto usando os potenciais secundários, os campos elétricos primários ( E^p) (Rijo, 2001) são calculados previamente e armazenados em arquivos numéricos. A determinação destes campos primários é realizada usando convolução numérica com os filtros de Anderson (Anderson, 1989).

A metodologia que este trabalho propõe é idêntica à metodologia proposta por Souza (Souza, 2007), uma vez que ambas aplicam o critério de Galerkin sobre as Equações (8) e (9). Porém, Souza (2007) aplica a identidade

apenas à esquerda da Equação (9), enquanto neste trabalho se aplica a identidade (Eq. (10)) tanto à esquerda quanto à direta da Equação (9), sendo esta a diferença fundamental entre ambas as metodologias. Ao se aplicar a identidade da Equação (10) àdireita da Equação (9), estamos admitindo a continuidade dadivergência de J^s.

A Equação (10) é uma identidade das funções bases (N_k e N_i) dos elementos finitos tridimensionais (nas direções y e z, a expressão é semelhante, utilizando derivadas parciais nas direções y e z, respectivamente) (Zienkiewicz & Taylor, 2000).No Apêndice A, as Equações (8) e (9) são solucionadas pelo método dos elementos finitos tridimensionais aplicando o critério de Galerkin, juntamente com a aplicação da identidade, Equação (10).

O problema inverso do mCSEM 3D

A detecção e/ou delineação de corpos resistivos tridimensionais usando como dados observados apenas a componente E_x do campo elétrico localizado no assoalho marinho, para apenasuma estação de fonte primária (DEH), não é simples. A literatura acerca do método mCSEM mostra coerentes resultados apenas quando são utilizadas, no processo de inversão 3D, as componentes x, y e z do campo elétrico total para várias estações de fonte primária (Zhdanov & Yoshioka, 2005). Desta forma, este trabalho também apresenta um algoritmo para inversão de dados eletromagnéticos 3D oriundos de fontes primárias galvânicas (DEH), usando apenas a componente E_x do campo elétrico (tanto a fase quanto a amplitude do referido campo) localizado no assoalho marinho, juntamente com a metodologia de modelagem direta do mCSEM 3D segundo os potenciais EM ( A,Φ) descrita anteriormente.

O problema da inversão geofísica aqui considerado consiste na minimização da função objetivo normalizada:

onde

desde que

Na Equação (11) p^s = é o vetor dos parâmetros geoelétricos parametrizados (Rijo & Almeida,2005), segundo com i = 1,2,3,¼,m (desde que > 0), na s-ésima iteração. O vetor das variáveis paramétricas na s-ésima iteração é dado por t^s = sendo m o número de parâmetros geoelétricos, n o número de observações com n > m e as componentes do vetor y^obs dos dados observados, tanto a amplitude quantoa fase, provenientes da componente E_x (in-line ) do campo elétrico medido no receptor localizado no assoalho marinho. As componentes f_k(p^s) são do vetor f(p^s), na s-ésima iteração, com k = 1,....,2n, tanto da fase quanto da amplitude dacomponente E_x do campo elétrico (in-line ) proveniente da modelagem direta do mCSEM 3D, também na s-ésima iteração.Cada ponto de medida do campo elétrico fornece sua amplitude e fase; logo para n pontos de medidas do referido campo tem-se n amplitudes e n fases, daí se ter 2n no primeiro somatório da Equação (11). O parâmetro l é o multiplicador de Lagrange e q é o número de parâmetros vinculados. O vetor γ = (γ_i(j)) com i ∈ {1,2,3,....,m} e j = 1,2,....,q contém os vínculos absolutos paramétricos. Por exemplo, digamos quesabemos os valores dos parâmetros p₂, p₈ e p_m os quais desejamos vinculá-los. Logo, q = 3, e assim através da parametrização acima sabemos que as variáveis paramétricas são t₂, t₈ e t_m, respectivamente, com γ₂₍₁₎ := t₂, γ₈₍₂₎ := t₈ e γ_m(3) := t_m. Por tanto, o vetor dos vínculos absolutosparamétricos é dado por

Finalmente o vetor v^s = com i ∈ {1,2,3,....,m} e j = 1,2,....,q contém as variáveis paramétricas vinculadas, na s-ésima iteração. Para o exemplo acima, o vetor das variáveis paramétricas vinculadas tem a forma

Sabendo que o vetor p^s está em função do vetor t^s e, pelo exposto acima, para cada conjunto específico de componentes discriminadas foi associado um vetor também específico. Considerando os vetores ω = (ω_k) e φ(t^s) = (Φ_k(p^s)), pode-se então expressar a Equação (11) de forma vetorial

No lado direito da igualdade da Equação (12) se tem a segunda parcela que representa a introdução de informação a priori . No presente caso esta informação é de caráter absoluto (Medeiros & Silva, 1996; Luiz, 1999).

O algoritmo de inversão que propomos é baseado no método de Marquardt (Press et al., 1992), juntamente com a estratégia proposta por Medeiros & Silva (1996) e Luiz (1999), que é exatamente a estratégia de introdução de informação a priori e que para o presente problema de inversão é especificamente chamada de introdução de vínculos absolutos. Desta forma, para a minimização da função objetivo, Equação (11) ou Equação (12), segundo o método de Marquardt (Press et al., 1992), tem-se a equação matricial iterativa

onde A_s é a matriz sensibilidade, cujos elementos são as derivadas parciais da função vetorial φ(t^s) em relação às variáveis paramétricas na aproximação t^s. Esta matriz possui a forma

e tem a ordem 2n×m. é a transposta da matriz A_s, λ é o parâmetro de Marquardt, e I é a matriz identidade de ordem m×m, M, de ordem m×m, é a matriz cujas componentes de suas linhas são formadas por zeros, exceto na sua componente localizada na diagonal principal, cujo índice é igual ao índice da variável paramétrica que se deseja vincular. Desta forma, esta componente possuirá o valor unitário. De acordo com o exemplo anterior sobre os vínculos absolutos, a matriz M apresentará a seguinte forma:

Assim, para a (s+1)-ésima iteração tem-se a aproximação do vetor das variáveis paramétricas dada por

No Apêndice B deste trabalho encontram-se uma abordagem teórica do método de Marquardt (Press et al., 1992) e o algoritmo de inversão de dados do mCSEM 3D, usando as equações apresentadas acima.

RESULTADOS

Resultados da modelagem direta do mCSEM 3D segundo os potenciais EM (A,Φ)

Para avaliarmos os sinais eletromagnéticos mCSEM 1D, 2D e 3D, faz-se necessário o uso de modelos geoelétricos em ambientes de prospecção de hidrocarbonetos em águas profundas.Desta forma consideremos, de modo comum a todos os modelos geoelétricos, um dipolo elétrico horizontal (DEH) 30 m acima do assoalho marinho, nas coordenadas x = 0 e y = 0, operando a uma frequência de 0,125 Hz. Inicialmente propomos um modelo geoelétrico de referência onde o ar possui uma resistividade elétrica infinita ρ₀ = ∞, o mar possui uma resistividade ρ₁ = 0,3 Ωm e espessura de h₁ = 1000 m e rochas sedimentares encaixantes de resistividade ρ₂ = 1 Ωm.

Em seguida, propomos um modelo geoelétrico com reservatório unidimensional, isto é, infinitamente longo nas direções x e y, em conformidade com o modelo de referência acima descrito, a menos de uma camada plano paralela representando um reservatório de hidrocarbonetos com resistividade ρ₃ = 100 Ωm, espessura de 100 m a uma profundidade em relação à superfície marinha h₂ = 1000 m (Fig. 3). Já para o modelo geoelétrico com reservatório bidimensional, propomos uma camada plano paralela de reservatório de hidrocarbonetos com resistividade ρ₃ = 100 Ωm, espessura de 100 m a uma profundidade em relação à superfície marinha h₂ = 1000 m e largura limitada na direção x (Fig. 1).

Finalmente para o modelo geoelétrico com reservatório tridimensional, propomos um volume tabular de reservatório dehidrocarbonetos com resistividade r₃ = 100 Ωm, espessura de 100 m a uma profundidade em relação ao assoalho marinhoh₂ = 1000 m, largura limitada na direção x e comprimento também limitado na direção y (Fig. 2).

Inicialmente comparamos as amplitudes da componente E_x (in-line ) provenientes dos modelos com hidrocarboneto 1D e 2D propostos, sendo tais amplitudes normalizadas pelas amplitudes das componentes E_x( in-line ) resultante do modelo de referência descrito no início desta seção.

O Experimento 2-2D (Fig. 4), corresponde ao modelo 2D acima caracterizado, com uma largura na direção x de 5 km. Já o Experimento 4-2D (Fig. 4), também corresponde ao modelo 2D acima caracterizado, porém com uma largura na direção xde 10 km. Os resultados 1D e 2D (Fig. 4) foram obtidos porRijo (2003, 2004), respectivamente.

Observa-se na Figura 4 a coerência física dos gráficos referentes ao Experimento 2-2D e Experimento 4-2D, uma vez que o primeiro advém de um resistor de apenas 5 km de largura e o segundo de um resistor de 10 km de largura, onde tais resistores propiciaram sinais de menor amplitude para o primeiro gráfico (Experimento 2-2D, curva: ¼) e de maior amplitude para osegundo gráfico (Experimento 4-2D, curva: ººº).

A curva tracejada (Fig. 4) possui duas funções: a primeira consiste em analisar a anomalia proveniente de um reservatório unidimensional, e a segunda, de servir como gráfico padrão e/ou de referência para avaliar propostas de soluções numéricas de modelagens do mCSEM 2D e 3D.

Agora compararemos as amplitudes da componente E_x ( in-line ) provenientes dos modelos 3D com hidrocarbonetos, sendo tais amplitudes normalizadas pelas amplitudes das componentes E_x (in-line ), geradas pelo modelo de referência anteriormente descrito.

As amplitudes normalizadas e observadas na Figura 5 provêm das simulações feitas através do algoritmo proposto neste trabalho, ou seja, segundo a metodologia de modelagem mCSEM 3D através dos potenciais EM (A,Φ) anteriormente descrita.

A simulação denominada de Experimento 2-3D (Fig. 5), corresponde ao modelo 3D acima caracterizado, com uma largura na direção x de 5 km e comprimento na direção y de 10 km. Já a simulação denominada de Experimento 4-3D (Fig. 5), também corresponde ao modelo 3D acima caracterizado, com uma largura na direção x de 10 km e comprimento na direção y de 10 km.

Finalmente comparamos as amplitudes da componente E_x (in-line ) provenientes dos modelos com hidrocarbonetos 1D, 2D e 3D propostos, sendo tais amplitudes normalizadas pelas amplitudes das componentes E_x ( in-line ) provenientes do modelo de referência. Mais especificamente, fazemos a comparação do campo normalizado do Experimento 2-2D (Fig. 5, curva: ....) em relação ao campo normalizado do Experimento 2-3D (Fig. 5, curva contínua) e também comparamos o campo normalizado proveniente do Experimento 4-2D (Fig. 5, curva: ººº) em relação ao campo normalizado proveniente do Experimento 4-3D (Fig. 5, curva: - · -).

Observa-se na Figura 5 as aproximações dos campos normalizados da componente E_x (in-line ) devido aos Experimentos 2-3D e 4-3D (curvas: contínua e - · -) do mCSEM 3D, em relação aos campos normalizados da componente E_x (in-line ) devido aos Experimentos 2-2D e 4-2D (curvas: ···· e ººº) do mCSEM 2D, respectivamente. Tanto as curvas segundo o mCSEM 2D, quanto do mCSEM 3D apresentam coerência com a variação da geometria do reservatório representado pelas variações de comprimento e largura. Observa-se a convergência das soluções do mCSEM 3D, propostas nestetrabalho, em relação às soluções do mCSEM 2D propostas por Rijo (2004).

Resultados da modelagem inversa do mCSEM 3D

Realizamos a inversão de dados do mCSEM 3D referente a um modelo específico, onde cada gama de dados do mCSEM 3D(os dados observados da componente E_x do campo elétrico total registrados nos receptores localizados no assoalho marinho) está relacionado a uma dada frequência. Desta forma, para realizarmos a inversão dos dados do mCSEM 3D, implementamos ao nosso código computacional, de inversão geofísica, as equações do problema direto do mCSEM 3D que foram modelados segundo os potencias EM (A,Φ) e solucionados segundo o método dos elementos finitos tridimensionais. Como parte complementar do código de inversão geofísica, implementamos a função objetivo normalizada (Eq. (11)), a ser minimizada segundo o método de Marquardt (Press et al., 1992).

O modelo geoelétrico tridimensional que iremos usar e queé mostrado na (Fig. 6), consiste em um corpo de hidrocarbonetos de resistividade elétrica ρ₃ igual a 50 Ωm, uma largura na direção x de 3 km, um comprimento da direção y também de3 km e com 50 m de espessura. Localizado a 1000 m de profundidade a partir do assoalho marinho e em meio a rochas sedimentares com resistividade elétrica ρ₂ igual a 1 Ωm. Com uma lâmina d'água igual a 1000 m de espessura com resistividade elétrica r₁ igual a 0,25 Ωm. O DEH (estático durante os dois processos de inversão) está localizado a 970 m da superfície do mar, nas coordenadas x = 0 e y = 0, com seu eixo nadireção x.

As seções do modelo geoelétrico acima, tanto através do plano z×x, quanto através do plano z×y são mostradas nas Figuras 7 e 8, respectivamente, a seguir.

As interposições de seções em profundidade do modelo geoelétrico (Fig. 6), através do plano x×y (horizontal), são mostradas na Figura 9 a seguir. Inicialmente, tem-se, a partir do fundo do mar: o dipolo elétrico horizontal (DEH) com seu eixo na direção x a uma altura do assoalho marinho de 30 m.Em seguida, para baixo, tem-se os receptores localizados e afixados no assoalho marinho e que estão representados pelo símbolo *, onde em cada linha de receptores, que são paralelas ao eixo x, os mesmos apresentam um espaçamento de500 m entre si. Desta forma a geometria de aquisição proposta por este trabalho é constituída por 170 receptores e um DEH localizado a 970 m da superfície do mar, nas coordenadas x = 0 e y = 0 (imóvel) e operando nas frequências de 0,075 Hz e 0,175 Hz. Finalmente tem-se mais uma seção horizontal que figura o corpo de hidrocarbonetos com uma profundidade de soterramento, em relação ao assoalho marinho, de 1000 m e em sua adjacência as rochas sedimentares encaixantes.

Apresentado o modelo geoelétrico a ser reconstituído através da metodologia de inversão de dados do mCSEM 3D proposta e descrita na segunda subseção da seção de metodologias, focalizaremos o referido modelo em apenas uma vizinhança do corpo de hidrocarbonetos (Fig. 10). Trata-se da região de interesse, onde tanto sua profundidade, quanto sua espessura são vinculadas, assim como células adjacentes de rochas sedimentares. Todos esses vínculos podem ser eventualmente implementados ao processo de inversão graças às interpretações de seções sísmicas em profundidade.

Agora, para realizar a primeira reconstituição e/ou inversão geofísica referente ao modelo da Figura 10, foi feita a simulação do modelo proposto (Fig. 6) com o DEH operando na frequência de 0,075 Hz, a fim de se possuir tanto as amplitudes quanto as fases da componente E_x do campo elétrico total capturado pelos receptores localizados no assoalho marinho. A estas amplitudes foram acrescidos 7% de ruído gaussiano e nas fases foram acrescidos 3,3% de ruído gaussiano. Esta etapa consiste na criação dos dados observados sintéticos a serem usados na primeira reconstituição.

Portanto, de posse dos dados observados sintéticos e aplicando a metodologia descrita na segunda subseção da seção de metodologias, tem-se a primeira reconstituição (Fig. 11), referente ao modelo geoelétrico caracterizado pela Figura 10.

Para este primeiro processo de reconstituição do referido modelo geoelétrico sob uma frequência de 0,075 Hz (Fig. 6), usamos o parâmetro de Marquardt igual a 10, o multiplicador de Lagrange igual a 1000, sendo que o tempo de duração do processo foi de 23 horas, com 107 iterações, e com o uso de um procedimento de parada que é apresentado na segunda subseção do Apêndice B deste trabalho. Este processo de parada consiste em atingir uma ordem menor que a tolerância de 10 através da norma do máximo sobre o vetor resultante da diferença entre dois vetores consecutivos dos parâmetros estimados. Tanto o parâmetro de Marquardt, quanto o multiplicador de Lagrange foram propostos de forma empírica, ou seja, para o parâmetro de Marquardt se admitiu o valor inicial igual a 10, já para o parâmetro de Lagrange foram admitidos vários valores positivos até quando se verificou tanto a convergência das estimativas dos parâmetros geoelétricos, quanto a minimização da função objetivo normalizada (Eq. (11)).

Finalmente, para realizar a segunda reconstituição referente ao modelo da Figura 10, foi feita a simulação do modelo proposto (Fig. 6), com o DEH operando na frequência de 0,175 Hz, a fim de se possuir tanto as amplitudes, quanto as fases, da componente E_x do campo elétrico total capturado pelos receptores localizados no assoalho marinho. A estas amplitudes foram acrescidos 7% de ruído gaussiano e nas fases foram acrescidos 3,3% de ruído gaussiano. Esta etapa consiste na criação dos dados observados sintéticos a serem usados na segunda reconstituição.

Portanto, de posse dos dados observados sintéticos e aplicando a metodologia descrita na segunda subseção da seção de metodologias, tem-se a segunda reconstituição (Fig. 12), referente ao modelo geoelétrico caracterizado pela Figura 10. Para esta segunda reconstituição do referido modelo geoelétrico sob uma frequência de 0,175 Hz (Fig. 6), usamos o parâmetro de Marquardt igual a 10, o multiplicador de Lagrange igual a 1000, sendo que o tempo de duração do processo foi de 49 horas para realizar 250 iterações sem o uso do procedimento de parada. Os parâmetros de Marquardt e de Lagrange foram obtidos de maneira idêntica ao primeiro processo de reconstituição descrito acima.

Como condição inicial para se dar início ao processo de reconstituição, onde os resultados finais são mostrados acima (Figs. 11 e 12), usou-se sobre cada célula do modelo inicial o valor da resistividade das rochas encaixantes igual a 1 Ωm, ou seja, a nossa condição inicial parte do pressuposto que nãoexiste um corpo de hidrocarbonetos. Porém, os resultados finais acima referenciados mostram a existência de um corpo de baixa condutividade elétrica, isto é, inferiores a 0,1 S/m. Para maiores frequências, a condição inicial acima descrita é inadequada, pois as curvas, tanto das amplitudes, quanto das fases dos dados observados, apresentam maiores distâncias entre as curvas das amplitudes e das fases, respectivamente, sendo estas últimas amplitudes e fases oriundas da condição inicial acima descrita,ou seja, de um modelo na ausência de hidrocarbonetos. Explicando de outra forma, o código de inversão geofísica que usa a derivada numérica, deve partir de uma condição inicial onde oseu modelo inicial calcula uma curva inicial ligeiramente próxima à curva dos dados observados, para assim viabilizar a minimização da função objetivo. Todas as duas reconstituições acima observadas foram calculadas no Laboratório PROEM/IG/ UFPA em uma máquina com CPU 3,00 GHz, 3,00 GB de RAMe administrada por um sistema operacional Linux. O primeiro processo de reconstituição, representado pela Figura 11, possui um critério de parada (Apêndice B), diferentemente do segundo processo, pois este realizou um total de iterações pré-estabelecidas de 250.

CONCLUSÕES

Entre os métodos geofísicos eletromagnéticos usados na exploração de hidrocarbonetos em águas profundas, exceto na zonado pré-sal, vem se destacando nos últimos dez anos o método mCSEM - marine Controlled Source Electromagnetic e que também vem recebendo grande atenção tanto pela indústria exploratória de petróleo quanto por centros de pesquisas. Devido à sua importância atual, faz-se necessária a sua pesquisa, tanto na prática quanto na simulação em laboratório de modelos geoelétricos de ambientes de hidrocarbonetos em águas profundas. Sendo assim, este trabalho contribuiu primeiramente tanto para a modelagem, quanto para a solução numérica do problema direto do mCSEM 3D. Em relação à modelagem, admitiu-se a continuidade da divergência do campo elétrico primário (Eq. (9)) de maneira que se tornou possível calcular a solução numérica do mCSEM 3D através do método dos elementos finitos tridimensionais segundo os potenciais EM (A,Φ), sem obter instabilidade numérica e incoerência na solução. Os experimentos desenvolvidos mostraram aproximações coerentes, no sentido qualitativo e quantitativo, quando avaliadas com as soluções numéricas devido aos experimentos propostos pela literaturado mCSEM 1D, 2D.

Este trabalho contribuiu também para inversão de dados do mCSEM 3D. Foi realizada a união, tanto de partes constituintes da primeira contribuição deste trabalho, quanto da segunda contribuição deste trabalho, a saber, a função objetivo normalizada proposta, Equação (11), a ser minimizada segundo o método de Marquardt (Press et al., 1992). Origina-se assim um código computacional que se propõe a realizar a inversão geofísica de dados do mCSEM 3D. Os resultados de inversão devido a este código computacional mostraram-se coerentes quando comparados ao modelo geoelétrico proposto (Fig. 6), de ambiente de hidrocarbonetos em águas profundas.

APÊNDICE A

APLICAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS TRIDIMENSIONAIS

Reescrevendo as Equação (8) e Equação (9) em termos de coordenadas cartesianas, obtemos o primeiro sistema de quatro equações a seguir Aplicando o critério de Galerkin (Becker et al., 1981) no primeiro sistema de quatro equações, Equações (A-1)-(A-4), e considerando as condutividades constantes em cada elemento, tem-se o segundo sistema de quatro equações a seguir

Aplicando o critério de Galerkin (Becker et al., 1981) no primeiro sistema de quatro equações, Equações (A-1)-(A-4), e considerando as condutividades constantes em cada elemento, tem-se o segundo sistema de quatro equações a seguir

em que N_k são funções bases, onde, para elementos tetraedrais, os valores de k variam de 1 a 4 (Zienkiewicz & Taylor, 2000).Usando a integração por partes (Teorema de Green), pode-se expandir as equações desse segundo sistema, Equações (A-5)-(A-8), formando o terceiro sistema de equações a seguir para k = 1,2,3,4.

para k = 1, 2, 3, 4.

Expandindo as componentes dos potenciais magnéticos e elétrico Φ^s juntamente com as componentes e do campo elétrico em termos das funções bases N_l para um tetraedro (Zienkiewicz & Taylor, 2000), onde os valores de l também variam de 1 a 4, tem-se

e

Se as expansões acima são substituídas nas equações do terceiro sistema acima, Equações (A-9)-(A-12), juntamente com as aplicações das condições de continuidade dos potenciais A^s e Φ^s tem-se o quarto sistema de equações a seguir

para k = 1,2,3,4.

A forma das equações do quarto sistema, Equações (A-13)-(A-16), resulta do cancelamento das integrais de superfície dos potenciais A^s e Φ^s, na formação da matriz global em suas respectivas contribuições internas.

O segundo membro, da última equação, Equação (A-16), do quarto sistema formado acima, Equações (A-13)-(A-16), pode ser expresso da seguinte forma:

Aplicando tanto no primeiro membro quanto no segundo membro, da última equação, Equação (A-16), do quarto sistema formado, as identidades

(Zienkiewicz & Taylor, 2000), tem-se para o segundo membro da última equação, Equação (A-16), do quarto sistema formado

Resulta que as integrais de superfície, nas identidades acima, também se cancelam na formação da matriz global em suas respectivas contribuições internas, devido ao fato de se considerar a continuidade da divergência do campo elétrico primário. Assim, a última equação, Equação (A-16), do quarto sistema, Equações (A-13)-(A-16), toma a forma

para k = 1,2,3,4. Essa expressão é diferente da usada por Souza (2007), dada por

para k = 1,2,3,4.

Para um elemento tetraedral genérico as funções bases são

com os índices k e l variando de 1 a 4 (Zienkiewicz & Taylor, 2000), onde n_e é o volume do tetraedro genérico e a_k, b_k, c_k, d_k,a_l, b_l, c_l e d_l são constantes a serem determinadas em termos das coordenadas x, y e z de cada vértice do tetraedro genérico. Finalmente, substituindo N_k e N_l nas equações do quarto sistema formado, Equações (A-13)-(A-16), obtém-se o seguinte sistemade equações lineares local para um elemento tetraedral genérico:

para k = 1,2,3,4.

Desta forma tem-se o sistema de equações lineares local, acima, composto por uma matriz de dimensão 16×16 e um vetor fonte de dimensão 16×1 para cada tetraedro pertencente ao domínio tridimensional da solução. Deste modo, cada sistema local é implementado computacionalmente para contribuir na formação do sistema de equações lineares global, sistema este associado a uma matriz complexa, simétrica e esparsa, que será resolvido numericamente pelo método do gradiente biconjugado de Souza et al. (2005).

APÊNDICE B

UMA ABORDAGEM TEÓRICA DO MÉTODO DE MARQUARDT

Sejam A ⊂ R^m um aberto e a aplicação U: → R de classe C². Consideremos o problema de minimizar a aplicação U(x) no aberto A. Desta forma, aproximemos U(x) pela série de Taylor até segunda ordem na vizinhança do ponto p^s∈ A, ou seja, estabelecemos a primeira aproximação a seguir: U( x)

em que

é o gradiente de U(x) no ponto p^s e

é a matriz Hessiana de U(x) no ponto p^s. Ao invés de H(p^s), consideremos a matriz Ha(p^s) como sua aproximação, que é constituída apenas de suas derivadas parciais de primeira ordem.

Agora, admitindo que o ponto p^s+1∈ A da vizinhança do ponto p^s é o minimante da aplicação U(x), isto é, ∇U(p^s+1) = 0, tomemos o gradiente da primeira aproximação acima contendo a substituição da matriz Hessiana pela matriz aproximada Ha(p^s):

E em seguida apliquemos, na segunda aproximação acima,o ponto p^s+1, resultando na terceira aproximação a seguir:

em que Δp = p^s+1 -p^s e Ha(p^s)^-1 é a matriz inversada matriz Hessiana aproximada.

Assim, o método de Levenberg-Marquardt consiste e madmitir:

1 ≤ i,j ≤ m, resultando na última aproximação a seguir:

Quando o parâmetro λ, chamado parâmetro de Marquardt, (λ > 0) é "grande" tem-se o domínio da matriz H'a(p^s) que contém sua diagonal dominante, admitindo assim matriz inversa, H'a(p^s)^-1, observado na última aproximação acima. Quando λ > 0 é "pequeno", tem-se o predomínio do gradiente ∇U(p^s), também observado na última aproximação acima.

Finalmente para minimizar U(x) no aberto A ⊂ R^m, usando a última aproximação acima, computemos inicialmente U(a^s), a^s∈ A, e tomemos um valor para λ ∈ [0.001,10]. Através da última aproximação acima calculemos Δa, fazemos a^s+1 := a^s + Δa e computemos U(a^s+1). Se U(a^s+1) ≥ U(a^s), então tomemos λ := λ*10 e através da última aproximação acima calculemos novamente Δa, tomemos a^s+1 := a^s + Δa e computemos novamente U(a^s+1).Se U(a^s+1) < U(a^s), então tomemos λ := λ/10 e a^s := a^s+1, voltando assim a computar U(a^s) e a calcular Δa atése atingir a tolerância estabelecida.

O ALGORITMO DO MÉTODO DE MARQUARDT COM VÍNCULOS ABSOLUTOS

Para minimizar a função objetivo normalizada U(p^s) (Eq. 11), utilizou-se o seguinte algoritmo listado abaixo:

ENTRADA

• As componentes com k = 1,¼,2n de dadosobservados (amplitude e fase) da componente E_x do campo elétrico total medido no receptor localizado no assoalho marinho;

• uma aproximação inicial dos parâmetros geoelétricos, do multiplicador de Lagrange l, do vetor γ = (γ_i(j)) com i Î {1,2,3,¼,m} ej = 1,2,....,q, dos vínculos absolutos paramétricos, do parâmetro de Marquardt λ, o número máximode iterações N e a tolerância TOL.

SAÍDA

Os parâmetros geoelétricos estimados , , ou uma mensagem de que o número máximo de iterações foi excedido.

Passo 1

Faça com i = 1,2,....,m. Observação: são as variáveis paramétricas. Faça chave := 1.

Passo 2

Faça s := 1.

Passo 3

Enquanto (s ≤ N) execute os Passos 4 a 11.

Passo 4

Se (chave = 1), então:

Passo 4.1

Calcule através do método direto do mCSEM 3D proposto as componentes f_k(p^s) com k = 1,2, 3,....,2n, formada tanto pelas amplitudes quanto pelas fases da componente E_x do campo elétrico total medidos no assoalho marinho. Observação: p^s = p^s-1.

Passo 4.2

Calcule Φ_k(p^s) = com k = 1,....,2n. Calcule U(p^s) e faça Ga := U(p^s).

Passo 4.3

Calcule (A_s)_(k,i) = com k = 1,....,2n e i = 1,2,....,m. Observação: t^s = t^s-1 e ts .

Passo 4.4

Calcule Δt^s (Eq. 13) e faça t^s+1 := t^s + Dt^s. Faça , com i = 1,2,....,m. Observação: .

Passo 5

Se (chave = 2), então:

Passo 5.1

Use os resultados calculados na iteração anterior dos Passos 4.1, 4.2 e 4.3 para se efetuar os cálculos do Passo 4.4, uma vez que se possui um novo valor para o parâmetro de Marquardt l.

Passo 6

Calcule através do método direto do mCSEM 3D proposto as componentes f_k(p^s+1), com k = 1,2, 3,....,2n, formada tanto pelas amplitudes quanto pelas fases da componente E_x do campo elétrico total medidos no assoalho marinho.

Passo 7

Calcule Φ_k(p^s+1) = f_k(p^s+1/y^obs_k, com k = 1,....,2n. Calcule U(p^s+1) e faça Gb := U(p^s+1). Faça erro := max|p^s+1 - p^s|.

Passo 8

Se (Gb < Ga), então:

Passo 8.1

Faça l := l*0.1, t^s₁ := t^s+1_i, com i = 1, 2,....,m, e com i = 1, 2,....,m. Faça chave := 1.

Passo 9

Se (Gb ≥ Ga), então:

Passo 9.1

Faça λ := λ*10, : = , com i = 1, 2,....,m, e , com i = 1, 2,....,m. Faça chave := 2.

Passo 10

Se (erro < TOL), então:

SAÍDA

('Parâmetros geoelétricos estimados',

SAÍDA

('Erro', erro). O procedimento foi bem-sucedido. PARE.

Passo 11

Faça s := s+1.

Passo 12

Se (s > N), então:

SAÍDA

('O número máximo de iterações foi excedido.'). O procedimento não foi bem-sucedido. PARE.

AGRADECIMENTOS

Os autores agradecem à ANP/PRH06 pela infraestrutura do Laboratório PROEM do IG/UFPA. O autor (F.L.A.) agradece ao CNPq pela bolsa de doutorado e, em particular, também agradece ao Victor C.T. Souza, da UFPa, pelas discussões esclarecedoras sobre o método dos elementos finitos tridimensionais. O autor (F.L.A.) também agradece in memoriam ao Prof. Luiz Rijo porter lhe orientado em sua tese de doutoramento.

Recebido em 19 abril, 2010 / Aceito em 15 julho, 2011

Received on April 19, 2010 / Accepted on July 15, 2011

NOTAS SOBRE OS AUTORES

Frayzer Lima de Almeida. Técnico em Mineração pela Escola Técnica Federal do Pará, Bacharel em Matemática, Mestre em Geofísica (Métodos Elétricos e Eletromagnéticos) e Doutor em Geofísica (Métodos Elétricos e Eletromagnéticos) (tese defendida e aceita em 05/02/2010), pela Universidade Federal do Pará. Atualmente é Prof. Adjunto da Faculdade de Matemática do Campus Universitário de Castanhal da Universidade Federal do Pará. Áreas de interesse: Otimização Matemática, Geofísica dos Métodos Elétricos e Eletromagnéticos e Estatística (Análise Multivariada e Inferência Bayesiana).

Luiz Rijoin memoriam . Graduado pela Universidade Federal de Pernambuco (UFPE) em Geologia, em 1965, e Matemática em 1969, PhD pela Universidade de Utah em 1977 com pós-doutorado no período de 1983/1984 pela mesma universidade, atuou como professor titular da Universidade Federal do Pará (UFPA), responsável pela disciplina de Métodos Geofísicos Eletromagnéticos no curso de pós-graduação em Geofísica, desde 1977 até o seu falecimento. Consultor ad hoc e membro de vários comitês do CNPq, Capes, Finep, Pronex, Adimb. Participou de diversas bancas examinadoras de concursos públicos e de defesa de teses e dissertações. Desenvolvia pesquisas na área dos métodos geofísicos eletromagnéticos aplicados à exploração de petróleo e gás. Foi durante muitos anos Pesquisador 1-A do CNPq, com mais de 60 trabalhos publicados em revistas científicas indexadas e em anais de congressos internacionais.

• ANDERSON WLO. 1989. A hybrid fast Hankel transform algorithm for electromagnetic modeling. Geophysics, 54: 263-266.
• BECKER EB, CAREY GF & ODEN JT. 1981. Finite elements - An Introduction. New Jersey: Prentice-Hall, 258 p.
• CONSTABLE SC, PARKER RL & CONSTABLE CG. 1987. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52: 289-300.
• CONSTABLE S & WEISS CJ. 2006. Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling.Geophysics, 71: 43-51.
• CONSTABLE S & SRNKA LJ. 2007. An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration.Geophysics, 72: WA3-WA11.
• CHRISTENSEN NB & DODDS K. 2007. Marine Controlled-Source Electromagnetic Methods 1D inversion and resolution analysis of marine CSEM data. Geophysics, 72: 27-38.
• LUIZ JG. 1999. Informação a priori na inversão de dados Magnetotelúricos. Tese de Doutorado, UFPA. 92 p.
• MacGREGOR L, SINHA M & CONSTABLE S. 2001. Electrical resistivity structure of the Valu Fa Ridge, Lau Basin, from marine controlled-source electromagnetic sounding. Geophys. J. Int., 146: 217-236.
• MEDEIROS WE & SILVA JBC. 1996. Geophysical inversion using approximate equality constrains. Geophysics, 61: 1678-1688.
• PRESS WH, TEUKOLSKY SA, VETTERLING WT & FLANNERY BP. 1992. Numerical Recipes in Fortran 77. Cambridge University Press. 2 ed., p. 678-680.
• RIJO L. 2001. Teoria dos Métodos Eletromagnéticos I, II e III. UFPa, Departamento de Geofísica. Disponível em: <http://www.rijo.pro.br> . Acesso em: 01 nov. 2006.
• RIJO L. 2003. Modelagem de dados MCSEM 1D. Relatório de Atividade do Convênio UFPA/PETROBRAS/FADESP.
• RIJO L. 2004. Modelagem de dados MCSEM 2.5D. Relatório de Atividade do Convênio UFPA/PETROBRAS/FADESP.
• RIJO L & ALMEIDA FL. 2005. Constrained 1-D inversion of MCSEM data on resistive oil reservoir. In: 9 International Congress of the Brazilian Geophysical Society, Salvador, Brazil. 1 CD-ROM.
• SOUZA VCT. 2007. Modelagem numérica de dados MCSEM 3D usando computação paralela. Tese de Doutorado, UFPA, 110 p.
• SOUZA VCT, RIJO L & SILVA MW C. 2005. The preconditioned biconjugate gradient algorithm applied to geophysical electromagnetic modeling. In: 9 International Congress of the Brazilian Geophysical Society, Salvador, Brazil. 1 CD-ROM.
• TARANTOLA A & VALETTE B. 1982. Generalized Nonlinear Inverse Problems Solved Using the Least Squares Criterion. Reviews of Geophysics and Space Physics, 20(2): 219-232.
• UEDA T & ZHDANOV MS. 2005. Fast numerical modeling of marine controlled-source electromagnetic data using quasi-linear approximation. In: SEG/Houston 2005 Annual Meeting, p. 506-509.
• ZHDANOV MS & YOSHIOKA K. 2005. Tree-dimensional iterative inversion of the marine controlled-source electromagnetic data. In: SEG/ Houston 2005 Annual Meeting, p. 526-529.
• ZIENKIEWICZ OC & TAYLOR RL. 2000. The finite element method -Volume 1: the basis. Oxford: Butterworth Heinemann, 707 p.

Datas de Publicação

• Publicação nesta coleção
  10 Fev 2012
• Data do Fascículo
  Jun 2011

Histórico

• Aceito
  15 Jul 2011
• Recebido
  19 Abr 2010