BRPI1106565A2 - Dispositivo e método para calcular convergências de ângulos de 3d a partir de migração de tempo reverso - Google Patents

Abstract

Dispositivo e método para calcular convergências de ângulos de 3d a partir de migração de tempo reverso. A presente invenção refere-se a um método para calcular convergências de imagens comuns de domínio de ângulo (adcig). O método inclui calcular uma fonte de campo de ondas pf de uma fonte sísmica: calcular um campo de ondas pb de receptor de receptor sísmico; aplicar um algoritmo de transformação de fourier antivazamento (alft) para transformar o campo de ondas de fonte pf para um domínio de número de ondas; aplicar o algoritmo de alft ao campo de ondas de receptor para transformar o campo de ondas de receptor no domínio do número de ondas; determinar uma condição de ondas de receptor no domínio do número de ondas; determinar uma condição de formação de imagem para a fonte de alft e campos de ondas de receptor no domínio de número de ondas; computar um ângulo de reflexão 0 e um ângulo de azimute 0 do campo de ondas de fonte pf e campo de ondas de receptor pb no domínio de número de ondas; calcular os adcig no domínio do número de ondas; e aplicar uma transformação de fourier rápida inversa (fft) para determinar os adcig no domínio de espaço

Description

Relatório Descritivo da Patente de Invenção para "DISPOSITIVO

E MÉTODO PARA CALCULAR CONVERGÊNCIAS DE ÂNGULOS DE 3D A PARTIR DE MIGRAÇÃO DE TEMPO REVERSO".

Remissão a Pedidos Relacionados O presente pedido reivindica a prioridade e benefício decorren- tes do pedido de patente provisório U.S. N°. 61/386.099, depositado em 24 de setembro de 2010, para "Convergências de Ângulos 3D a Partir de Mi- gração de Tempo Reverso", cujo inteiro teor fica incorporado na sua totali- dade neste contexto por referência.

Antecedentes Campo Técnico As concretizações da matéria em apreço exposta neste contexto referem-se de uma maneira geral aos métodos e sistemas e, mais em parti- cular, aos mecanismos e técnicas para proporcionarem de forma mais efici- ente uma imagem precisa de uma estrutura do subsolo.

Discussão Sobre os Antecedentes Durante os ultimo anos, o interesse no desenvolvimento de no- vos campos de produção de petróleo e gás aumentou drasticamente. Entre- tanto, a disponibilidade de campos de produção terrestres é limitada. Desta forma, a indústria estendeu agora a perfuração a localizações em alto-mar que parecem conter uma vasta quantidade de combustível fóssil. A perfura- ção em alto-mar é um processo dispendioso. Assim, aqueles que se dedi- cam à perfuração em alto-mar necessitam de saber onde perfurar a fim de evitar um poço seco. A aquisição e processamento de dados sísmicos marítimos ge- ram um perfil (imagem) da estrutura geofísica sob o fundo do mar. Muito embora este perfil não proporcione uma localização exata para o petróleo e gás, ele sugere, para aqueles treinados no campo, a presença ou ausência de petróleo e/ou de gás. Desta forma, proporcionar-se uma imagem de alta resolução das estruturas sob o fundo do mar constitui um processo em an- damento.

Durante um processo de convergência sísmica, tal como ilustra- do na Figura 1, uma embarcação 10 arrasta uma série de detectores acústi- cos proporcionados em arrastes 12. Os arrastes podem ser dispostos hori- zontalmente, ou seja, dispostos a uma profundidade z1 constante em rela- ção a uma superfície 14 do oceano. A embarcação 10 também arrasta um conjunto de fonte de som 16 que é configurado para gerar uma onda acústi- ca 18. A onda acústica 18 propaga-se descendentemente no sentido do fundo do oceano 20 e penetra no fundo do oceano até eventualmente uma estrutura refletora 22 (refletor) refletir a onda acústica. A onda acústica refle- tida 24 propaga-se ascendentemente até ser detectada por meio do detector 12. Os dados registrados são então processados para produzirem uma i- magem exata da subsuperfície. O processamento inclui várias fases, por exemplo, determinação de modelo de velocidade, pré-empilhamento, migra- ção, pós-empilhamento, e outros, que são conhecidos na técnica e, deste modo, a sua descrição é omitida neste contexto. O progresso na formação de imagem de profundidade de pré- empilhamento foi considerável no passado. O progresso teórico proporcio- nou melhores métodos para extrapolação de campos ondulatórios medidos na superfície da terra para a subsuperfície, e o progresso prático foi vincula- do às migrações mais estreitamente com a construção e interpretação do modelo de velocidade. A migração é o processo de propagação, por exem- plo, um campo ondulatório medido em um local de receptor para um refletor localizado na subsuperfície. A migração também pode ser aplicada aos campos ondulatórios gerados por meio de uma origem.

Nas áreas de subsuperfícies complexas, as dificuldades de for- mação de imagens são devidas a dois componentes: construção do modelo de velocidade de profundidade de pré-empilhamento e algoritmos de migra- ção. A construção do modelo de velocidade avalia o modelo de velocidade (por exemplo, como a onda sonora se propaga através das várias camadas da terra) para a simulação da propagação de ondas sísmicas que ocorrem durante a migração. Este modelo forma a parte de comprimento de onda longa (macro) do modelo de terra, e a migração proporciona a parte de com- primento de onda curta (refletividade). A tomografia baseada em raios sís- micos é uma ferramenta amplamente usada para a construção de modelos, mas a não linearidade e incerteza dos algoritmos da tomografia baseada em raios expõe a tomografia como um elo fraco no processo de formação de imagem.

Outro elo fraco no processo de formação de imagem é a ilumi- nação sísmica deficiente das regiões por baixo do terreno detrítico complexo (por exemplo, sais, estruturas de arrastamento, e assemelhados), que difi- cultam ou tornam mesmo impossível a formação de imagem adequada. A iluminação deficiente do terreno detrítico é causada freqüentemente pela aquisição sísmica inadequada, por exemplo, por aquisição de arraste de a- zimute estreito 3D convencional quando é necessária aquisição de azimute amplo 3D.

Até recentemente, a migração de Kirchhoff foi o método principal usado para a migração de profundidade de pré-empilhamento. Este método demonstrou ser bem sucedido ao longo de numerosos exemplos quando as variações de velocidade são mínimas. Este método também formou a base para a migração de "amplitude real". Este algoritmo faz migrar os dados sísmicos de entrada um rastreio de cada vez ou um grupo local de rastreios de cada vez; estes processos significam que o custo de migração de Kirc- hhoff é proporcional ao número de rastreios de entrada. Não obstante, quando o número de rastreios de entrada é relativamente pequeno dentro da abertura de migração (como é usualmente o caso com os levantamentos de azimutes estreitos marinhos), a migração de Kirchhoff produz um algoritmo eficiente. Pelo contrário, quando o número de rastreios de entrada é grande dentro da abertura de migração (tal como é usualmente o caso com os le- vantamentos de azimutes amplos marinhos), a eficiência pode ser perdida quando a tarefa computacional se torna exigente.

Por outro lado, a utilização do traçado de raios para aproximar a função de Green da propagação de ondas pode comprometer a exatidão da migração de Kirchhoff, especialmente quando o campo ondulatório é com- plexo. Uma aproximação tradicional, por exemplo, a escolha da chegada de um único raio do campo ondulatório complexo em cada local da imagem, de- termina uma imagem com ruído em áreas onde existem muitas chegadas de raios. Os algoritmos de migração de Kirchhoff de múltiplas chegadas supe- ram este problema, mas eles tendem a ser complicados e relativamente ine- ficientes em 3 dimensões (3D).

De acordo com outras abordagens, os métodos de migração de profundidade por pré-empilhamento de feixes aproximam a função de Green com uma expressão que permite que várias chegadas de imagens sejam re- presentadas sem complicação algorítmica excessiva, e podem ser aplicados em um sentido de amplitude real. Não obstante, da mesma forma que ocor- re para a migração de Kirchhoff, a aproximação da função de Green usada por meio de migração de feixes baseia-se no rastreamento de raios e pode tornar-se imprecisa se o modelo de velocidade de migração contiver varia- ções extremamente fortes (por exemplo, corpos de sal) e precisar de suavi- zação excessiva. Ainda assim, a capacidade da migração de feixes repre- sentar estruturas complexas e de controlar determinados tipos de ruído de migração pode usualmente assegurar imagens significativamente melhores nas áreas complexas do que os algoritmos de migração de Kirchhoff de che- gada simples.

Muito embora os métodos de migração de Kirchhoff e de feixe u- tilizem raios para proporcionar a aproximação das funções de Green para propagação de ondas, os chamados algoritmos de migração da equação de ondas utilizam as funções de Green de forma de onda plena que são gera- das numericamente, por exemplo, por meio de diferenças finitas. Os algo- ritmos computacionalmente mais eficientes para que isto seja feito são cha- mados coletivamente de migração de equação de direção única (OWEM).

Estes algoritmos decompõem os campos de ondas sísmicos dentro da terra em ondas ascendentes e ondas descendentes na suposição de nenhuma interação entre estes dois campos de ondas; ou seja, nenhuma geração de onda de retorno e de reflexão vertical durante a síntese da propagação de ondas. Sobre um corpo muito grande e crescente de exemplos, a OWEM solucionou os problemas de múltiplas chegadas melhor do que a migração de Kirchhoff de chegada única. Para os levantamentos sísmicos de azimute estendido, em que o número de rastreios de entrada é amplo em compara- ção com a abertura de migração, a OWEM tende a ganhar eficiência em re- lação à migração de Kirchhoff. Para esses levantamentos, a implementação eficiente dos algoritmos de OWEM pode ser construída seja para migração de tiro comum ou para migração de ondas planas. Não obstante, existem algumas limitações importantes dos algo- ritmos de OWEM. Primeira, as ondas de viragem estão ausentes na síntese de propagação de ondas, o que resulta em os eventos de alta inclinação em torno de 90° serem insuficientemente representados; segunda, a síntese de propagação de ondas somente assegura a exatidão da fase do campo de ondas, enquanto que as amplitudes do campo de ondas são muito menos confiáveis e precisam de correção adicional. O uso da equação de onda de duas vias na migração em pro- fundidade começou algum tempo atrás em um algoritmo chamado de migra- ção de tempo inverso (RTM). Não obstante, esta abordagem foi limitada de- vido à sua necessidade de alimentação de computador. Com os aumentos na energia de computador, a RTM foi rapidamente desenvolvida durante os últimos anos, e vantagens teóricas tais como propagação de ondas precisas de inclinação ilimitada e amplitudes aperfeiçoadas proporcionaram, na práti- ca, benefícios de formação de imagens. Por exemplo, nas áreas de flanco de pré-sal e de sal, as funções numéricas de Green provenientes da diferen- ça finita da equação de ondas de duas vias têm melhor comportamento de amplitude, de maneira que é mais fácil incorporar correções de amplitude na RTM do que na OWEM. Adicionalmente à sua capacidade de manipular dis- tribuições de velocidades complexas, muitos algoritmos de RTM atuais po- dem manipular meios anisotrópicos, tais como isotropia transversal vertical (VTI) e isotropia transversa inclinada (TTI). Em exemplos de representação de imagens de dados reais, a RTM de TTI proporcionou as melhores ima- gens em um exame topográfico de azimute amplo complexo do Golfo do México, embora os modelos de velocidade para a migração de TTI fossem simplificados como isotropia transversalmente conformável estruturalmente (STI), que requer o eixo simétrico anisotrópico compatível com os vetores normais dos refletores.

Com a exatidão aperfeiçoada da RTM vem a sensibilidade au- mentada para a exatidão do modelo de velocidade. Esta sensibilidade oca- siona melhoramento notável nas imagens de RTM quando o modelo de ve- locidade é preciso, mas ela também provoca degradação digna de nota das imagens de RTM quando o modelo de velocidade não é preciso. Por esta razão, a análise da velocidade de migração é mais importante para a RTM do que ocorre para outros métodos de migração de profundidade. A ligação entre migração e construção de modelo de velocidade é um conjunto de convergências de imagens comuns (CIGs) produzidas por meio dos algoritmos de migração. Uma CIG é um conjunto de imagens, to- das na mesma localização de imagem (usualmente na localização do refletor na subsuperfície), com cada imagem formada a partir de diferentes subcon- juntos de dados de entrada. Por exemplo, um volume único de dados de deslocamento comum / azimute comum, que é um subconjunto do conjunto de dados sísmicos de pré-empilhamento adquirido pleno, pode ser usado para reproduzir a terra em 3D. A coleta de imagens a partir de todos os subconjuntos de dados com diferentes formas de deslocamento e azimute proporciona a formação das CIG. As CIG incluem múltiplos rasteios. As CIG podem ter todos os rastreios com diferentes deslocamentos (com a tota- lidade da informação azimutal somada), ou as CIG podem ser dotadas de todos os rasteios com diferentes deslocamentos e azimutes.

As CIG são usadas comumente para análise de variação de am- plitude de domínio de profundidade com deslocamento (AVO), e análise de velocidade baseada em migração. Com um modelo de velocidade correta, todas as imagens na mesma localização de imagem devem-se concentrar sob a mesma profundidade, fazendo com que os eventos de reflexão nas CIG pareçam planos. O achatamento dos eventos sísmicos nas CIG é um dos critérios para a validação do modelo de velocidade por meio de análise de focalização. Quando os eventos nas CIG não são planos, os geofísicos aperfeiçoam os seus modelos de velocidade de migração por meio de análi- se da curvatura dos eventos, utilizando a análise para orientar a atualização de velocidade.

Para a migração de Kirchhoff, não há custo adicional significativo para computar CIGs de deslocamento comuns (COCIGs). Por outro lado, a migração de volumes comuns - deslocamento por meio de OWEM ou RTM é dispendiosa, de maneira que as COCIG não estão de uma maneira geral disponíveis para esses métodos de migração.

As limitações de qualidade das COCIG são causadas, em parte, pelas limitações subjacentes da migração baseada em raios. De maneira mais fundamental, as COCIG sofrem dos artefatos de migração devido aos múltiplos caminhos de propagação das ondas, sejam os métodos de migra- ção capazes ou não de manipular múltiplos caminhos de campos de ondas com exatidão, provocando potencialmente dificuldades para a análise de ve- locidade e análise de amplitude contra ângulo de reflexão (AVA). Com efei- to, as CIG cujos rasteios são indexados por meio de qualquer atributo na su- perfície de registro, tais como deslocamento de fonte / receptor ou ângulo de incidência de superfície da energia de origem, são suscetíveis a tais artefa- tos. Sob este aspecto foi demonstrado que uma condição necessária para CIG isentas de artefatos é a de serem parametrizadas em um domínio de ângulo de subsuperfície, tais como em um ângulo de reflexão ou ângulo de abertura. Isto foi ilustrado em 2D utilizando-se migração de Kirchhoff de múltipla chegada no conjunto de dados sintéticos de Marmousi e trabalho subseqüente ampliou isto mostrando meios anisotrópicos, ou 3D utilizando- se as CIG no ângulo de reflexão / ângulo de azimute, e para análise em 3D em domínio de múltiplos ângulos (ângulo de reflexão, ângulo de inclinação, ângulo de azimute, e assemelhados).

Comparadas com as migrações de feixe e de Kirchhoff de múlti- plas chegadas, as OWEM e RTM parecem ter capacidades limitadas para as CIG indexadas no domínio de deslocamento de superfície. No domínio de ângulo de superfície, propôs-se uma abordagem que emite como saída CIGs de deslocamento de subsuperfície locais a partir de OWEM e então converte as mesmas em CIGs de domínio de ângulo (reflexão) de subsuperfície (AD- CIGs). A conversão das CIGs de deslocamento de subsuperfície locais em ADCIGs tem uma forma simples no caso isotrópico de 2D. Esta abordagem requer que a condição de representação de migração seja aplicada em uma faixa de deslocamentos de subsuperfícies, formando as CIG de deslocamen- to de subsuperfície; em seguida, uma transformação de Fourier 2D é aplica- da à CIG de deslocamento local; então o número de ondas de transformação é mapeado para o ângulo de reflexão. O procedimento é realizado conver- gência por convergência localmente, e é um algoritmo razoavelmente efici- ente. Não obstante, a fórmula de conversão de convergência é com- plexa em 3D uma vez que ela requer a produção de CIGs de deslocamento de subsuperfície locais indexadas por meio de deslocamento X e desloca- mento Y, e uma transformação de Fourier 5D para as CIGs de deslocamento de subsuperfície locais para o domínio de número de ondas, seguida por um mapeamento de alto dimensional. Para a implementação, mesmo produzin- do-se um vetor de deslocamento de subsuperfície local (X, Y) as CIGs a par- tir da OWEM são computacionalmente intensivas. Uma abordagem alternati- va consiste em produzir as CIGs de retardamento de tempo, e então conver- te-las para ADCIGs. A CIG de retardamento de tempo tem apenas um eixo adicional nas imagens (tempo de retardamento). Uma saída da CIG de re- tardamento de tempo é 4D (x, y, z, tempo de retardamento). Não obstante, para dados sísmicos de azimute estendidos, as ADCIG 3D esperadas são elas próprias 5D (x, y, z, e ângulo incidente / de reflexão, ângulo de azimu- te), de modo que a condição de formação de imagem de deslocamento de tempo não proporciona informação suficiente para uma conversão plena em ADCIGs de 3D quando aplicada a dados sísmicos de azimute estendidos.

Como um resultado, de uma maneira geral, às convergências de ângulos finais falta informação azimutal. Além disso, a conversão de um domínio de dados para outro sofre de problemas de amostragem, degradando prova- velmente a resolução das CIG finais.

As abordagens expostas anteriormente para se obterem as AD- CIG para OWEM e RTM são indiretas, prosseguindo através de um domínio intermediário (deslocamento de subsuperfície ou atraso de tempo de ima- gem). É igualmente possível computar as ADCIG por meio de decomposi- ção direta dos campos de ondas na subsuperfície em seus componentes di- recionais locais. Assim sendo, este trabalho foi realizado em algoritmos de migração isotrópicos em 2D ou 3D, e produz ADCIGs somente no ângulo de reflexão (não azimute). Além disso, esta abordagem foi aplicada somente aos propagadores de campos de ondas de uma via, tornando difícil reter in- formações de amplitude confiáveis. Muito embora haja outras diferenças en- tre os propagadores para OWEM e RTM, não existe diferença estratégica essencial nos algoritmos para a saída de CIG.

Por essa razão, seria desejável proporcionarem-se sistemas e métodos que evitassem os problemas e inconvenientes descritos anterior- mente.

Sumário De acordo com uma concretização exemplificativa, proporciona- se um método para calcular, em um dispositivo de computação, convergên- cias de imagens comuns de domínio de ângulo (ADCIG) para o uso na gera- ção de uma imagem final imagem de uma subsuperfície da terra. O método inclui uma etapa de calcular com propagação antecipada de campo de on- das, em um domínio de espaço de tempo, uma fonte de campo de ondas pF de uma fonte sísmica; uma etapa de calcular com propagação de retorno de campo de ondas, no domínio de espaço de tempo, um campo de ondas pB de receptor de um receptor sísmico configurado para detectar uma onda sísmica gerada pela fonte sísmica; uma etapa de aplicar um algoritmo de uma transformação de Fourier anti-fugas (ALFT) para transformar o campo de ondas de fonte pF para um domínio de número de ondas; uma etapa de aplicar o ALFT algoritmo ao campo de ondas de receptor para transformar o campo de ondas de receptor no domínio de número de ondas; uma etapa de determinar uma condição de formação de imagem para a fonte de ALFT e campos de ondas de receptor no domínio do número de ondas; uma etapa de computação de um ângulo de reflexão Θ e um ângulo de azimute φ do campo de ondas de fonte pF e campo de ondas de receptor pB no domínio de número de ondas; uma etapa de calcular as ADCIG no domínio de núme- ro de ondas; e uma etapa de aplicar uma transformação de Fourier rápida inversa (FFT) para determinar as ADCIG no domínio de espaço.

De acordo com outra concretização exemplificativa, proporciona- se um meio capaz de ser lido por computador que inclui instruções executá- veis por meio de computedor, em que as instruções, quando executadas, implementam um método para calcular convergências de imagens comuns de domínio de ângulo (ADCIG) para o uso na geração de uma imagem final de uma subsuperfície da terra. O método inclui as etapas observadas ante- rio rmente.

De acordo ainda com outra concretização exemplificativa, pro- porciona-se um dispositivo de processamento configurado para calcular con- vergências de imagens comuns de domínio de ângulo para o uso na geração de uma imagem final de uma subsuperfície da terra. O dispositivo de pro- cessamento inclui uma interface configurada para receber to receber dados sísmicos e um processador conectado à interface. O processador é configu- rado para calcular com propagação antecipada de campo de ondas, em um domínio de espaço de tempo, um campo de ondas pf de fonte uma fonte sísmica, calcular com propagação de retorno de campo de ondas, no domí- nio de espaço de tempo, um campo de ondas pB de receptor de um receptor sísmico configurado para detectar uma onda sísmica gerada pela fonte sís- mica, aplicar um algoritmo de transformação de Fourier anti-vazamento (ALFT) para transformar o campo de ondas de fonte pF para um domínio de número de ondas, aplicar o algoritmo de ALFT ao campo de ondas de recep- tor para transformar o campo de ondas de receptor no domínio de número de ondas, determinar uma condição de formação de imagem para a fonte de ALFT e campos de ondas de receptor no domínio de número de ondas, computar um ângulo de reflexão Θ e um ângulo de azimute φ do campo de ondas de fonte pF e campo de ondas de receptor pB no domínio de número de ondas, calcular as ADCIG no domínio de número de ondas, e aplicar uma transformação de Fourier rápida inversa (FFT) para determinar as ADCIG no domínio de espaço.

Descrição Breve dos Desenhos Os desenhos anexos, que são incorporados no relatório e consti- tuem uma parte do mesmo ilustram uma ou mais concretizações e, em con- junto com a descrição, expõem estas concretizações. Nos desenhos: A Figura 1 é um diagrama esquemático de uma configuração de aquisição de dados convencional. A Figura 2 é um diagrama esquemático que indica a definição de um ângulo de reflexão Θ e um ângulo de azimute cp de subsuperfície local. A Figura 3 é um fluxograma que ilustra um método para calcular convergências de imagens comuns. A Figura 4 é um fluxograma que ilustra uma transformação de Fourier anti-vazamento. A Figura 5 é um gráfico que mostra uma função de teste. A Figura 6 é um gráfico de uma função que reconstitui a função de teste da Figura 5 por meio de FFT. A Figura 7 ilustra espectros obtidos para a função da Figura 5 por meio de FFT e ALFT de acordo com uma concretização exemplificativa. A Figura 8 ilustra os espectros da Figura 7 em maiores detalhes de acordo com uma concretização exemplificativa. A Figura 9 ilustra a função da Figura 5 reconstruída de acordo com uma concretização exemplificativa. A Figura 10 mostra a função reconstruída e a função original e o seu agrupamento errôneo de acordo com uma concretização exemplificativa. A Figura 11 ilustra a localização dos receptores ande uma fonte de som de acordo com uma concretização exemplificativa. A Figura 12 ilustra dados sísmicos sintéticos de acordo com uma concretização exemplificativa.

As Figuras 13 e 14 ilustram localizações de imagens de acordo com uma concretização exemplificativa. A Figura 15 ilustra as CIG de ângulo de azimute comum calcula- das com base nos dados sintéticos da Figura 12 de acordo com uma concre- tização exemplificativa. A Figura 16 ilustra uma pilha de todas as contribuições a partir de uma aquisição simulada de azimute ampla completa de acordo com uma concretização exemplificativa. A Figura 17 ilustra uma imagem empilhada de uma linha central de acordo com uma concretização exemplificativa. A Figura 18 ilustra ADCIGs em 2D convertidas a partir das con- vergências de deslocamento de subsuperfície de acordo com uma concreti- zação exemplificativa. A Figura 19 ilustra ADCIGs 3D de azimute zero de acordo com uma concretização exemplificativa. A Figura 20 mostra uma seção empilhada em linha TTI RTM em uma área complexa. A Figura 21 mostra uma única ADCIG 3D a partir de uma migra- ção de um subconjunto de azimute estreito dos dados de azimute estendido de acordo com uma concretização exemplificativa. A Figura 22 mostra uma CIG de ângulo 3D a partir de uma mi- gração de um conjunto de dados de azimute estendido pleno de acordo com uma concretização exemplificativa. A Figura 23 é um fluxograma de um método para calcular con- vergências de ângulo 3D a partir de migração de tempo reverso de acordo com uma concretização exemplificativa; e A Figura 24 é um diagrama esquemático de um dispositivo de processamento para calcular convergências de ângulo 3D a partir de migra- ção de tempo reverso de acordo com uma concretização exemplificativa.

Descrição Detalhada A descrição seguinte das concretizações exemplificativas refere- se aos desenhos anexos. Os mesmos números de referência em desenhos diferentes identificam os mesmos elementos ou elementos assemelhados. A descrição detalhada seguinte não é limitativa da invenção. Em vez disso, o escopo da invenção é definido pelas reivindicações em anexo. As concreti- zações seguintes são discutidas, por razões de simplicidade, com relação à terminologia do campo de ondas de fonte, campos de ondas de receptor, convergências de imagens em comum, e migração de tempo reverso para o processamento de dados sísmicos. Não obstante, as concretizações que serão discutidas em seguida não ficam limitadas a estes métodos e concei- tos, mas podem ser aplicadas a outros métodos para produzir imagens finais da subsuperfície. A referência por todo o relatório a "uma concretização" ou "a concretização" significa que um aspecto, estrutura, ou característica particu- lar descrita em conexão com uma concretização está incluído em pelo me- nos uma concretização da matéria em apreço exposta. Desta forma, o apa- recimento das frases "em uma concretização" ou "numa concretização" em vários locais por todo o relatório não está necessariamente referindo-se à mesma concretização. Além disso, os aspectos, estruturas ou característi- cas particulares poderão ser combinados de qualquer uma maneira adequa- da em uma ou mais concretizações.

De acordo com uma concretização exemplificativa, a RTM pro- duz CIGs isentas de artefatos tanto para atualização de velocidade quanto análise de amplitude, e estas CIG deverão conter, não suprimir, informações azimutais. Essas CIG pode ser indexadas por meio de ângulo de reflexão t ângulo de azimute, e as CIG pode ser obtidas por meio da realização de uma decomposição direcional de campo de ondas em localizações na sub- superficie. Estas CIG, chamadas CIG de domínio de ângulo 3D, podem ser produzidas em um meio isotrõpico ou em um meio anisotrópico STI.

De acordo com outra concretização exemplificativa, os campos de ondas de fonte propagadas para a frente e campos de ondas de receptor propagadas para trás são transformadas para o domínio do número de on- das. Deste modo, os vetores de números de ondas espaciais proporcionam a direção de fase da propagação de ondas. De acordo com uma aplicação, com a finalidade de reduzir o custo numérico, as CIG são formadas em uma janela local, equivalente a uma decomposição de ondas de plano local.

Ainda de acordo com outra concretização exemplificativa, utilize- se uma transformação de Fourier anti-vazamento (ALFT) para substituir a transformação de Fourier rápida (FFT). O algoritmo ALFT ajuda a reduzir o fenômeno de Gibbs que ocasiona sérios problemas de exatidão quando ope- ra em uma janela local pequena. Além disso, a ALFT delineia as ondas de plano local por meio da sua energia, que podem ser usadas para limitar o número de ondas de plano mais energético na condição de formação de i- magem de convolução, em vez de usar o conjunto completo das mesmas.

Este conjunto ajuda a reduzir de forma notável o custo de computação, tal como discutido em seguida.

Antes de se discutirem os novos métodos, acredita-se que seja apropriada uma revisão da formulação para a migração de amplitude real de 3D com saída de convergência de ângulos de subsuperfície. Esta formula- ção proporciona uma formula assintótica que é usada posteriormente para aplicação à RTM. Para uma equação de onda acústica, a refletividade R na localização x de subsuperfície (os símbolos em negrito representam um ve- tor) depende de um ângulo de reflexão Θ e de um ângulo de azimute cp de subsuperfície local tal como ilustrado na Figura 2. Desta forma, a refletivida- de R depende da localização x, do ângulo de reflexão e do ângulo de azímu- te como R = R(x, θ, φ). A Figura 2 também mostra uma posição s de uma fonte e uma posição r de um receptor. A fonte em s pode ser qualquer fonte conhecida, por exemplo, fonte de pistola ou fonte de EM, e o receptor em r pode ser qualquer receptor conhecido, por exemplo, geofone, hidrofone, e assemelhado. Os vetores de números de ondas ks e kr são definidos como as direções de fase de propagação de ondas na localização de imagem x para os campos de ondas provenientes da localização de fonte s e localiza- ção de receptor r respectivamente. Os pares de ângulos (0S, cps), (ΘΓ, cpr), e (9m, <Pm) são, respectivamente, os ângulos de inclinação e os ângulos de a- zimute dos vetores ks, kr, e a sua soma k de vetores, respectivamente. Por meio da lei de Snell, o plano 30 ortogonal ao vetor k é um plano refletor po- tencial (que coincide com um refletor real quando ks e kr formam um par de números de ondas especulares). As relações dos ângulos são definidas como: (1) O vetor unitário é definido como nx = (1, 0, 0). Na equação (1), as definições dos ângulos Θ e φ são para um caso geral, que é aplicável aos meios VTI e TTI isotrópicos. A aproximação de Kirchhoff para a síntese sísmica na modelagem para a frente estabelece que, no domínio de fre- qüência, o campo de ondas registrado õG(r, ω, s) na localização r devido a uma fonte na localização s, é dado por: <5G(r,ã),s) = Jí/xR(x, -(v Ga (r, ω, x)C?Ô (x,®,s) — G0 (r, ω, x)V G„ (x,6>,s)) em que G0(x, ω, s) é a função de Green no domínio de frequên- cia para propagação no meio de baixa prioridade a partir da localização s para a localização x, a operação matemática indica o conjugado comple- xo, R(x, θ, ψ) = R(x, θ, φ)η com R(x, θ, φ) sendo a refletividade dependente de ângulo e sendo η o vetor normal ao plano de reflexão, e sendo V é o ope- rador do gradiente espacial. A integral sobre o volume especial pleno é a soma das contribuições de refletividade diferencial proveniente de todas as localizações de subsuperfície x. Assintoticamente, pode-se expressar a fun- ção de Green como G0(r,a>,x) = A(r,x)eímT{T-*\ Gg(x,6),s) = Ã(x,s)e^lmr(x,’\ ^ e o gradiente da função de Green assintótica é dada por VG0(x,a>, s) = ~ia>V T(x, s)A(x, s)e~ía,T(x’’) ^ Nas equações (3) e (4), T(x, y) e A{x, y) são o tempo de percur- so e amplitude de raio a partir da localização y até x, respectivamente. Co- mo um resultado, pode-se expressar o campo de ondas refletidas na equa- ção (2) como: em que o núcleo integral de Kirchhoff K(r, x, s, ω) = iu>nq(r, x, s)A(r, x, s) contém um termo de vetor de inclinação de migração q(r,x,s) = VT(r,x) + V7\x,s). No caso isotrópico, também se pode expressar o vetor de inclinação de migração como q(r,x,s) - —n, em que Θ ê o v(x) ângulo de incidência, e v(x) é a velocidade de propagação de onda na loca- lização de imagem x. Este vetor de inclinação é combinado com um termo de produto de amplitude z((r,x,s) = ^(r,x)^(x,s) no núcleo de migração.

Observa-se que a equação (5) é válida somente com a aproximação assintó- tica de alta freqüência. O termo de fase é a soma dos tempos de percurso dos campos de ondas no ponto de imagem a partir da localização de fonte e da localização de receptor: 2"(r,x,s) =r(r,x) + T(x,s) (g) Se for observado um sismógrafo <5Gobl,(r,a>,s) em uma localiza- ção r devido a uma fonte na localização s, estes dados observados deverão ser compatíveis com os dados de reflexão sintetizados <SG(r5®,s). Esta coe- rência requer que o modelo de velocidade de baixa prioridade e o modelo de refletividade sejam os dois corretos. Esta suposição de migração de ampli- tude verdadeira implica que a velocidade de baixa prioridade obtida propor- cione uma estimativa razoável para as variações de comprimento de onda longas do modelo terra, de maneira que o tempo de percurso da propagação de onda é bem definido. A tarefa da inversão é encontrar um bom modelo de refletividade para se conseguir um ajuste entre os dados SGobs(r,a>,s) ob- servados e os dados <S<G(r,ft>,s) sintéticos. Utilizando-se um esquema de inversão de mínimos quadrados ponderados, pode ser gerada uma função objetiva para medir um desajuste entre os dados calculados õG{r,co,s) e os dados adquiridos/observados SGobs(r,co,s) como se segue: (7) em que |]*|j indica um modelo 1_2 e Q = Q($.,0,s,r) é o peso para a inversão dos mínimos quadrados. Este peso precisa ser determinado. A redução ao mínimo da função de desajuste (7) proporciona a solução dos mínimos quadrados para o modelo de refletividade Κ(χ,θ,φ). Esta solução é a chamada imagem de migração de amplitude real ou preservada, ou ima- gem de migração/inversão.

Partindo-se da equação (5), observa-se que os dados sintéticos õG(r,a>,s) dependem linearmente do modelo de refletividade Κ(χ,θ,φ), e a equação (7) coloca o problema a ser solucionado em um enquadramento de teoria de inversão linear. A resolução da equação (7) quanto à refletividade Κ{*,θ,φ) sob ângulos θα e φΰ específicos, produz R(x,6]0,φ0)· O produto R(x,90,<p0) pode ser realizado iterativamente utilizando-se um método Gauss-Newton, cuja primeira solução de iteração (supondo-se uma estimati- va de refletividade inicial nula) proporciona o modelo de refletividade depen- dente de ângulo: (8) No lado direito da equação (8), o primeiro termo é o inverso do Hesseno e o Segundo termo é o gradiente da função objetiva. O gradiente da função objetiva tem a expressão: (9) em que as duas funções S (δ(θ-θ0) e δ(φ-φα)) são funções de armazenamento de ângulos que restringem a integração para o ângulo de reflexão θ0 desejado e o ângulo de azimute φα. O Hesseno tem a 0YnrpsRãn· (10) em que x é a localização de foco real e Xo é a localização de fo- co com o modelo de velocidade de baixa prioridade. A segunda etapa da equação (10) é o resultado da mudança das variáveis de integração partin- do-se das coordena^00 oi,r'erfície (s, r) para os ângulos de subsuperfície (θ, φ) e a Jacobiana I o resultado desta mudança de variáveis. A aproximação na equação (10) é válida quando a aproximação assintótica de alta freqüência ω(Τ(r, x, s) —Τ(χ, x0, s)J« k ■ (x - x0) se mantém.

Quando o modelo de velocidade de baixa prioridade fica suficientemente próximo do modelo de velocidade verdadeira, as funções de Green em x0 e x ficam muito próximas dentro da banda de freqüências sísmicas. Portanto, na aproximação de primeira ordem sob uma aproximação assintótica de alta freqüência, o núcleo de amplitude poderá ser aproximado como: 3C(r,x,s,o) « fJC(r,x0,s,<») . (11) Então, enquanto o peso dos mínimos quadrados satisfizer (12) o lado direito da equação (10) pode ser proporcionado assintoti- camente equivalente à função delta de Dirac: (13) Neste caso, a Hessiana torna-se uma matriz diagonal: Η(χ,χ0ΑΛ) = <5·(χ-χ0). (14) A Jacobiana na equação (10) inclui a determinante de Beylkin.

Esta determinante proporciona a conexão entre os vetores de direção locais em uma localização de imagem e a geometria de aquisição real. Para uma aquisição de azimute ampla, nos casos de meios isotrópicos, ela pode ser desenvolvida como: em que v(x) é a velocidade no ponto de imagem x e p\ é o componente vertical do vetor de lentidão da propagação de onda na superfí- cie. Com a ajuda da equação (15), o peso Q da inversão dos mínimos qua- drados torna-se: (16) Observe-se que esta função ponderai não depende diretamente da amplitude do campo de ondulações. Esta produção regular provém da geometria de aquisição de azimute ampla. De acordo com um aspecto ge- ral, o peso depende, através da determinante de Beylkin, do desvio da geo- metria de aquisição em relação ao caso de deslocamento estendido, de azi- mute estendido ideal. A aquisição de azimute estendido cobre (em princípio) a superfície da terra completamente com um arranjo de receptores de área para registrar a energia proveniente de cada disparo e um arranjo de dispa- ros de área para enviar a energia a ser registrada em cada receptor. Como um resultado disto, o peso (16) perdeu a dependência na amplitude do cam- po de ondas que é usado para compensar a aquisição insuficiente, de ma- neira que é diferente do peso de migração de deslocamento comum tradi- cional. A substituição das equações (9), (14), e (16) na equação (8) proporciona a seguinte forma assintótica para migração de amplitude real de azimute estendido 3D com as ADCIG de saída: (17) que proporciona migração de Kirchhoff de amplitude real de 3D para ADCIGs tanto no ângulo de reflexão quanto no ângulo de azimute. Não obstante, a formula de inversão acima tem uma singularidade em θ = 0°, que é diferente da sua contraparte de inversão em 2D.

Para compreender porque razão acontece esta singularidade, supõe-se que os ângulos de reflexão e de azimute são distribuídos em uma semi-esfera como latitude e longitude, respectivamente. Θ = 90” dá o equa- dor, com todos os azimutes distribuídos igualmente no mesmo. θ = 0° é o pólo norte em que o ângulo de azimute é definido de forma ambígua. O ter- mo l/sin<9 na condição de formação de imagem supra é usado para com- pensar a diminuição da área quando o ângulo de reflexão tende para 0°. Não obstante, isto causa dificuldades na obtenção de amplitudes precisas em torno do ângulo de reflexão de 0o. A razão desta diferença se deve à função de armazenamento direto δ(φ-φ0).

Determina-se agora a migração de amplitude real de 3D para migração de tempo reversa. A expressão de migração de amplitude real (17) proporciona uma receita para produzir ADCIGs de 3D no ângulo de re- flexão e no ângulo de azimute. Não obstante, esta expressão é uma ex- pressão assintótica, e ela pode ser aplicada diretamente apenas aos algo- ritmos de migração que calculam explicitamente termos de amplitude e ba- seia-se em direções iniciais especificadas, ou seja, em migrações de Kirc- hhoff e de feixe. Para se aplicar a teoria à RTM, é necessário traduzir a ex- pressão (17) para uma formula de RTM de amplitude real. A fórmula de RTM resultante será ainda assintótica em termos de amplitude, mesmo que o motor de propagação de onda subjacente não o seja. Para simplificar a sua derivação, considera-se a RTM como isotrópica (o caso anisotrópico po- de ser tratado de forma assemelhada, com complicações extras na notação). É igualmente considerado decompor-se a migração plena na expressão (17) em migração tiro por tiro.

Se as múltiplas reflexões são ignoradas, o operador de migração de RTM é linear para os dados de entrada SGobs, e a saída de azimute es- tendida final é a soma das contribuições provenientes de todas as experiên- cias de tiro. Para produzir uma fórmula de amplitude real para RTM, as ex- pressões assintóticas para a função de Green (termos de amplitude e tempo de percurso) são substituídas pelas suas contrapartes exatas. Isto é conse- guido ao reescrever-se a expressão de migração (17) como: (18) Neste caso, R (x, θ0, φα) = R(x, θ0, φ0 ) /|qj = R(x, θ0, φ0 )ν(χ) / 2 cos θ„ é introduzido para manter uma consistência com uma definição de saída em um método tradicional. A integral sobre as localizações da fonte assegura que todos os rasteios sísmicos registrados contribuam para a imagem final e para as convergências de ângulos. Os campos de ondas pF e pB são da- dos por meio de: (19) A propagação para a frente da assinatura de fonte (campo de ondas de fonte), pF satisfaz: (20) Na equação (20), a integração de tempo da ondulação de fonte /(/) incorpora o termo pl em pF. A retro propagação (campo de ondas de receptor) dos dados registrados pB satisfaz à equação: (21) Nas equações (20) e (21), v = v(x) é a velocidade, S(x-s)f(t) é a assinatura de fonte, e v2 é o operador de Laplacian. A integral de equação final (18) sobre as localizações de tiros s é a soma de todas as contribuições para todos os tiros. Para um único tiro, a contribuição para as CIGs de ângulo 3D final é (22) Das duas condições de formação de imagem de migração co- nhecidas, correlação cruzada (multiplicação dos campos de ondas) e des- convolução (divisão dos campos de ondas), acredita-se que a condição de formação de imagem de correlação cruzada é apropriada para as ADCIG de 3D. A equação (22), e a sua contraparte teórica de raios (17), têm as duas uma forma de uma condição de formação de imagem de correlação cruzada, que é fácil de aplicar e numericamente estável.

Na expressão (22), os campos de ondas propagados para a frente e os campos de ondas propagados para trás estão em uma forma ge- ral, aplicável a qualquer migração baseada em equação de ondas, OWEM ou RTM. A fórmula de migração tem uma forma simples que é independente no modelo de velocidade, isotrópico ou anisotrópico STl.

Para se criarem as ADCIG a partir de RTM, a equação (18) ne- cessita dos campos de ondas de subsuperfície pFepB a partir das equa- ções (20) e (21). Uma vez que estes campos de ondas sejam computados, as direções de propagação de ondas (locais) ks e kr são determinadas em cada localização de imagem. Então a equação (1) proporciona o ângulo de reflexão Θ e o ângulo de azimute φ. A condição de formação de imagem na freqüência / domínio do número de ondas para o empilhamento pleno é dada por meio de: (23) Esta equação expressa a refletividade para cada ângulo de a- bertura / azimute como uma soma das contribuições de campos de ondas sobre uma faixa de ângulos de inclinação / azimute que obedece às leis de reflexão. As contribuições provenientes dos ângulos de inclinação / azimute que não concordam com essas contribuições de um refletor real tendem a se anular pela soma integrante. Com a finalidade de aperfeiçoar a eficiência dos cálculos, esta expressão é reescrita no domínio de número de ondas como: A soma sobre a freqüência na equação (24) proporciona uma condição de formação de imagem apropriada para a RTM. A equação (25) é usada para computar a refletividade dependente de ângulo por meio de uma transformação de Fourier inversa especial de 3D.

As etapas seguintes ilustradas na Figura 3 resumem o novo pro- cedimento discutido anteriormente. Na etapa 300, calculam-se os campos de ondas pF e pB utilizando-se um motor de propagação de ondas de RTM.

Os campos de ondas provenientes da localização de fonte são uma propa- gação para a frente (em tempo) de uma condição de fonte de limite; os cam- pos de ondas provenientes das localizações de receptores são uma propa- gação para trás (em tempo) dos dados sísmicos adquiridos pré-filtrados.

Todos os campos de ondas podem ser armazenados em um disco local. Na etapa 302, uma transformação de Fourier espacial / temporal 4D é aplicada aos campos de ondas tanto de fonte quanto de receptor (para a frente em tempo para o campo de ondas de fonte, para trás em tempo para o campo de ondas de receptor) para expressar os campos de ondas no domínio de freqüência / número de ondas. Esta transformação de 4D pode incluir várias outras etapas que são discutidas mais adiante. Na etapa 304, a condição de formação de imagem (24) é aplicada aos campos de ondas para decompor a imagem por meio do ângulo de reflexão Θ e ângulo de azimute φ, de acordo com a equação (1). Na etapa 306, uma transformação de Fourier espacial inversa de 3D é aplicada a lt(u.,0o,<po) para obter ADCIG de 3D Rs(x,0o,tpo) utilizando-se a equação (25). O algoritmo proposto é computacionalmente intenso porque a condição de formação de imagem (24) é uma convolução sobre os vetores de ondas tridimensionais provenientes tanto do lado de fonte quanto do lado de receptor. Dentro do laço de convolução, proporciona-se um processo de armazenamento com armazenamento de ângulo de azimute e ângulo de re- flexão. Este armazenamento requer a computação dos ângulos de reflexão e de azimute θ0β φϋ para cada par ks, kr no laço interno. Com o lado ex- terno de frequências, o número total de laços computacionais para a forma- ção de imagem de azimute plena pode ser sete para uma migração de tiro comum. O esforço computacional para a propagação de ondas para este procedimento de convergência de ângulo é o mesmo que para a RTM con- vencional. Não obstante, as transformações de Fourier de 4D e a decompo- sição de ângulo no domínio de número de ondas adicionais aumentam o custo computacional porque a condição de formação de imagem é uma con- volução multidimensional. Para uma RTM de tiro comum que efetua uma formação de imagem de azimute largo no Golfo do México, as aberturas de migração são tipicamente de 20 km nas direções tanto em linha, quanto em linha cruzada, com uma taxa de amostragem em linha de 25 m e uma taxa de amostragem de linha cruzada de 40 m, e a profundidade máxima é tipi- camente de 15 km, com uma taxa de amostragem de profundidade de 15 m. O número total de números de ondas para um único campo de ondas de subsuperfície pode alcançar 4xl08 (krorks = 800x500x1000), de modo que a convolução de dois campos de ondas pode alcançar l.óxlO17 operações de pontos flutuantes. Além disso, um laço externo adicional está presente para calcular a imagem final e este laço externo pode ter aproximadamente 200 freqüências. Desta forma, o número total de operações de pontos flutu- antes pode ser 3.2 χ 101S. Para essa grande número de operações de pontos flutuantes, o algoritmo descrito anteriormente pode ir bem além da capaci- dade dos atuais supercomputadores.

De acordo com uma concretização exemplificativa, para se re- duzir o custo de computação, poderá ser aplicado um par de outros algorit- mos e/ou técnicas tais como discutidos em seguida. Por exemplo, relações de dispersão são usadas nos vetores de ondas para restringir a decomposi- ção e convolução. Para modelos de velocidade isotrópicos homogêneos, as normas de todos os vetores de lentidão são um valor constante fixo, permi- tindo a redução de sete laços para cinco. Para modelos de velocidade hete- rogêneos e/ou anisotrópicos, as normas dos vetores de lentidão não são mais um valor fixo quando a velocidade varia na janela local, mas elas po- dem ser restringidas até faixas situadas entre valores mínimos e máximos, determinados pelos valores máximo e mínimo da velocidade de fase. Nos casos com uma ampla faixa de velocidades (tais como ambientes marinhos sobre sal), a redução de custo computacional devido a esta restrição é muito menos dramático do que nos casos com uma pequena faixa de velocidades.

Então, os campos de ondas da imagem são decompostos em janelas espaciais locais sobrepostas. Nesta decomposição, tipicamente as dimensões de janelas locais têm uma dimensão igual nas três dimensões espaciais. Uma dimensão de janela de uns poucos comprimentos de ondas pode ser usada sob uma freqüência dominante com 50% de sobreposição em cada dimensão espacial. A utilização de janelas locais reduz o custo da convolução de 3D em várias ordens de grandeza uma vez que ele reduz o número total de números de ondas. A utilização de janelas locais também permite um uso mais eficiente das relações de dispersão porque o limite lo- cal das normas de vetor de lentidão é usualmente muito menor do que o limi- te global. Por exemplo, pela seleção de kr or ky = 64 χ64 χ64 e levando-se em consideração o número de frequências situar-se em torno de 200 e o número de janelas ser situado em torno de 1664 vezes 8, em que 8 é o nú- mero de janelas de sobreposição, o número total de operações de pontos flutuantes pode ser reduzido para cerca de 2,0 χ 1017. Não obstante, duas questões se colocam em decorrência do es- quema de janelas locais discutido anteriormente. Primeira, se a dimensão de janela é muito pequena, o pequeno número de amostragens espaciais na janela de entrada resulta em um pequeno número de números de ondas no domínio de Fourier. A representação esparsa das direções de propagação de ondas causada pelo pequeno número de números de ondas causará, por sua vez, resolução deficiente no domínio de ângulos. Por exemplo, 32 nú- meros de ondas na direção x resultarão em uma resolução de ângulo maior do que 6o.

Segunda, quando a dimensão da janela se torna menor, os efei- tos de truncamento de limite da janela tornam-se maiores, de forma que os artefatos de Gibbs podem contaminar os componentes do número de ondas e degradar os resultados da migração. Estes dois fenômenos foram investi- gados pelos inventores da presente invenção e foi proposto um esquema de sobre amostragem de domínio de número de ondas dentro de uma transfor- mação de Fourier anti-vazamento (ALFT) para aliviar estas duas questões. A ALFT com esquema de número de ondas de sobre amostra- gem proporciona boa resolução angular e reduz os artefatos de contorno de Gibbs. O algoritmo de ALFT pode superar o "vazamento espectral" quando os coeficientes de Fourier não são ortogonais depois da sobre amostragem de números de ondas, e é capaz de manusear dados desalinhados. Sucin- tamente, o algoritmo de ALFT opera como se segue. Primeiro, utiliza trans- formações de Fourier distintas (DFT) ou o seu algoritmo alternativo para promover uma estimativa de coeficientes de Fourier para todos os números de ondas k. Então, ele seleciona aquele com a energia máxima e adiciona o componente de energia máxima ao número de ondas de saída. Em seguida, ele realiza uma transformação de Fourier distinta reversa (DFT) de volta à grade de entrada. Finalmente, o algoritmo de ALFT subtrai os dados trans- formados de volta à etapa anterior a partir da entrada na grade de entrada.

Repetem-se estas etapas para cada um dos coeficientes de Fourier até que todos os coeficientes sejam avaliados.

Desta forma, o algoritmo de ALFT pode ser implementado nas etapas seguintes tal como se encontra ilustrado na Figura 4. Na etapa 400, inicializam-se todos os componentes de Fourier da DFT para zero. Na etapa 402, computam-se todos os componentes de Fourier dos dados de entrada.

Observa-se que poderão ser usadas outras transformações de Fourier em vez da DFT. Na etapa 404, seleciona-se o coeficiente de Fourier com a e- nergia máxima, e acumula-se a sua contribuição para o componente de Fou- rier existente. Na etapa 406, subtrai-se a contribuição deste coeficiente a partir dos dados de entrada. Então, na etapa 408 determina-se se todos os coeficientes de Fourier foram calculados. Se a resposta é sim, então pára- se a ALFT de outro modo repetem-se as etapas 402 a 406 até os dados de entrada residuais atualizados serem suficientemente pequenos.

Tal como observado anteriormente, de acordo com uma concre- tização exemplificativa, um esquema de sobre amostragem de domínio de número de ondas é aplicado dentro do algoritmo de ALFT. De acordo com uma aplicação, uma relação de sobre amostragem de 8 pode ser usada em todas as três dimensões espaciais. O esquema de sobre amostragem au- menta o número de números de ondas drasticamente. Depois de sobre a- mostragem dos números de ondas para cada campo de ondas em uma jane- la local, o número de números de ondas poderá ser tão alto quanto 224, rein- troduzindo assim a possibilidade de um alto custo de convolução. Por e- xemplo, tomando-se uma relação de sobre amostragem de 4, isso significa que krorky= 64x64x64x4 com um número total de operações de ponto flu- tuante de cerca de S.OxlO20, que constitui novamente um aumento no custo da computação.

Recapitulando-se, o uso da ALFT permite controlar o custo da convolução. O esquema de sobre amostragem na ALFT favorece o uso de pequena dimensão de janela local sobre a preservação de toda a informação de números de ondas a partir dos campos de ondas, e a pequena dimensão de janela tende a proporcionar uma melhor aproximação de ondas planas para os campos de ondas por meio de um esquema de sobre amostragem dentro do algoritmo ALFT. Além disso, a ALFT pode ser usada para revelar os números de ondas mais energéticas.

Em outras palavras, em aplicações reais, um número relativa- mente pequenos de números de ondas contribuirá com a maior parte da e- nergia do campo de ondas em uma janela local. Desta forma, uns poucos números de ondas energéticas provenientes tanto da fonte quanto dos cam- pos de ondas de receptor contribuirão para a convolução. Por exemplo, o número de ondas mais energéticas pode ser limitado a 50, ou seja, kr or k5 = 50 que produz um número total de operações de ponto flutuante de cerca de 6,656x109, que é um número de operações muito reduzido que po- de ser resolvido pelos computadores dos dias de hoje.

Para exemplificar a potência do novo método, uma função sim- ples y = O.bc traçada na Figura 5 é usada para ilustrar a solução da ALFT. A

Figura 5 expõe uma amostragem regular com localizações de amostragem xe [0,---,63]. Os dados amostrados incluem 64 amostras com uma taxa de amostragem de 1, ou seja, uma grade regular. Como os dados de entrada são amostrados em uma grade regular, os componentes de Fourier são or- togonais uns aos outros e não se proporciona qualquer vazamento de fre- qüência. Uma ALFT de avanço direto não pode aperfeiçoar a reconstrução neste caso ideal. Não obstante, com a taxa de amostragem do número de ondas da teoria FFT, a Figura 6 mostra a reconstrução continua dos compo- nentes de Fourier computados por meio da FFT. Os artefatos de toque que aparecem nas duas bordas 50a e 50b são o conhecido fenômeno de Gibbs.

Este exemplo indica que a truncagem da faixa de amostragem introduz ruído perto do limite, especialmente para aqueles campos de ondas cujos núme- ros de ondas não estão na localização exata de amostragem de número de ondas da teoria de FFT.

De acordo com uma concretização exemplificativa, os dados no domínio de número de ondas são sobre amostrados tais como discutidos an- teriormente, uma vez que na ALFT o número de ondas poderá ser um valor arbitrário. Em uma aplicação de janela de pequena dimensão, o número de rastreios de entrada para ALFT é menor do que para aplicações de janela de grande dimensão, especialmente no caso de numerosas dimensões. Esta nova abordagem economiza o custo computacional da ALFT. Não obstante, depois da sobre amostragem do número de ondas, o número dos números de ondas pode ser maior do que quando se utiliza a teoria de FFT. O núme- ro dos números de ondas pode inclusive alcançar o mesmo número dos nú- meros de ondas que quando se utilizam aplicações de janelas de grandes dimensões. Desta forma, esta abordagem resulta em alto custo computa- cional para a convolução de domínio do número de ondas. Não obstante, utilizando-se dimensões de janelas menores, os eventos sísmicos ficam mais próximos dos eventos lineares e deste modo, o subconjunto de "número de ondas pontiagudas" do número de ondas total poderá representar bem os dados de entrada. Esta vantagem está ilustrada por um exemplo de 1D nas Figuras 7 e 8. A Figura 8 mostra uma vista mais detalhada do espectro para um baixo número de ondas k, por exemplo, me- nor do que 0,1. As Figuras 7 e 8 mostram uma linha 60 correspondente a um espectro computado por FFT. Outro espectro computado por ALFT está ilustrado pela linha 62. Observa-se que a linha 62 descreve melhor a fre- quência esperada para o modelo. No exemplo da Figura 8, mostra-se um fenômeno interessante na área do primeiro número de onda não zero. O espectro de ALFT mostra somente um valor 64 com uma grande energia en- quanto a FFT ilustra claramente o vazamento de frequência. A Figura 9 mostra a função continua reconstruída (linha 70) da Figura 5. Comparada com a solução exata (linha 72) ilustrada na Figura 10, a função reconstruída 70 tem pequenos erros e as maiores discrepâncias nas duas áreas de limite 74 e 76 são como ilustradas na Figura 10.

Nas áreas complexas, o número de eventos pontiagudos no do- mínio de número de ondas é de uma maneira geral maior do que um. Desta forma, na prática, pode ser suficiente tomar apenas 50~100 dos números de ondas de energia mais alta em uma janela local durante uma fração de tem- po fixa em vez dos 224 números de ondas por convolução. Este esquema pode reduzir o custo de convolução de forma significativa, e também pode reduzir o custo de aplicar a condição de formação de imagem expressa pela equação (24). Na prática, o esforço computacional de exibição parcial, apli- cação da ALFT, e realização da formação de imagem de ângulo pleno, ex- cede notavelmente o esforço de uma RTM convencional (5~10 vezes mais dispêndio computacional), mas permanece viável com os atuais super com- putadores (por exemplo, agrupamentos de PC).

Os novos métodos discutidos anteriormente são agora aplicados a situações da vida real. Primeiro, a RTM em 3D é testada em um meio de velocidade constante. É utilizado um modelo que é dotado de cinco refleto- res planos com um espaçamento de profundidade de 1,0 km. Um conjunto de dados de tiro comum de azimute amplo é sintetizado por meio da simula- ção de um único registro de tiro de azimute amplo indicado pelo número de referência 92 na Figura 11. As localizações dos receptores são distribuídas em uma grade de 100 m por 100 m com um deslocamento máximo de 5,0 km nas duas direções X e Y tal como ilustrado também na Figura 11. Os dados foram sintetizados com aquisição de um único tiro localizado na posi- ção de superfície (5,0 km, 5,0 km), que é o ponto central do arranjo receptor 94.

Cinco telhas 100a-e de dados sísmicos sintéticos estão ilustra- das na Figura 12. Realizou-se uma RTM de tiro comum para emitir como saída CIGs de ângulo de 3D com oito ângulos de azimute. Selecionaram-se duas localizações de imagem (3,0 km, 3,0 km) e (3,75 km, 5,0 km) para emi- tir como saída ADCIGs de azimute pleno. As Figuras 13 e 14 ilustram a exa- tidão do modelo supra e a 15 mostra que os términos teóricos dos eventos são compatíveis com os resultados. Estas figuras também ilustram a resolu- ção do novo algoritmo. Uma vez que esta é uma RTM de um único tiro co- mum, as contribuições para as convergências de ângulos são um ponto para cada localização dentro do azimute correto. Estas são chegadas especula- res identificadas sob os ângulos de reflexão corretos pelo novo procedimen- to. A nova abordagem de computação de ADCIG de 3D ilustra uma resolução superior, a energia especular está bem concentrada para a posi- ção de ângulo teórico, e isto é bem resolvido tanto para o ângulo de reflexão quanto para o ângulo de azimute. Em seguida, é discutida com referência à Figura 16 uma aquisição de azimute amplo completo simulado por meio da translação lateral do registro de tiro único sobre uma grade regular de locali- zações de tiros, espaçados em uma grade de 100 m por 100 m regular sobre todo o levantamento. A Figura 16 mostra o empilhamento de todas as con- tribuições provenientes desta aquisição de azimute amplo completo. Os términos das convergências de ângulos sob as diferentes profundidades são compatíveis com a solução teórica, com boa preservação de amplitude.

Para formação de imagem estrutural complexa, algoritmos de ADCIG para OWEM ou RTM consideraram até agora somente ângulos de reflexão de 2D ou 3D, sem considerar o azimute de reflexão em 3D. De a- cordo com uma concretização exemplificativa, proporciona-se um exemplo de CIGs de domínio de ângulo de reflexão / azimute a partir de RTM utili- zando-se o conjunto de dados sintéticos de Sigsbee. O modelo original é 2D. Isto foi generalizado para um modelo de 2,5D pelo preenchimento do modelo na direção Y. Com a finalidade de manter válidas as amplitudes de 3D, utilizou-se um esquema de modelagem de FD de alta ordem para gerar dados sintéticos de 3D com aquisição de azimute zero (arraste único). Para este conjunto de dados de 3D restringidos, espera-se que as ADCIG de 3D tenham notavelmente energia com ângulos de reflexão plenos somente no azimute zero. A Figura 17 mostra a imagem empilhada da linha central (por exemplo, uma linha de controle) utilizando-se a RTM de 3D. Selecionaram- se dezoito localizações (representadas por meio de 18 linhas 110 na Figura 17) para emitir as CIG como saídas. As ADCIG 2D que foram convertidas a partir das convergências de deslocamentos de subsuperfície estão ilustradas na Figura 18. O deslocamento de subsuperfície máximo é de 5000 pés com um incremento de deslocamento de 75 pés. As ADCIG são amostradas a cada 2 graus. Estas convergências exigem um efeito de borrado. Este efei- to é devido a questões de amostragem na conversão de convergências do espaço para o número de ondas para o ângulo. As ADCIGs de 3D de azi- mute zero estão ilustradas na Figura 19. As CIG de ângulo 3D exibem me- lhor resolução de sinal para ruído (S/N) e resolução mais alta com termina- ções de eventos mais nítidas sob mais ângulos do que os resultados de 2D decorrentes da abordagem convencional.

De acordo com outra concretização exemplificativa, realizou-se um teste de azimute largo real a partir de uma pesquisa em águas profundas na área dos Garden Banks no Golfo do México. Os dados sísmicos foram adquiridos com deslocamento máximo em linha de 8,0 km e deslocamento de linha transversal de 4,0 km. A Figura 20 mostra uma seção empilhada em linha de RTM TTI em uma área complexa área que inclui uma mini bacia entre dois corpos de sal complexos com um topo irregular de geometrias de sal e saliências. Os eventos do sal básico estão bem definidos. Realizaram- se duas RTM de TTI.

Primeiro, migrou-se um subconjunto de azimute estreito dos da- dos de azimute largo. A Figura 21 mostra uma ADCIG 3D única proveniente desta migração. A faixa de ângulo de reflexão é de 0-40° com uma taxa de amostragem de 2o. Utilizaram-se seis setores de azimute de superfície com uma taxa de amostragem de 30°. As contribuições dos ângulos de azimute para as CIG também foram traçadas acima da ADCIG 3D. Para o subcon- junto de azimute estreito dos dados, a energia na ADCIG é concentrada no azimute zero, com um certo vazamento para os ângulos de azimute próxi- mos sob profundidades rasas. Segundo, fez-se migrar o conjunto de dados de azimute estendido pleno. A Figura 22 mostra o CIG de ângulo 3D. A dis- tribuição de energia é distribuída mais uniformemente sobre os ângulos de azimute, exceto sob grandes profundidades em torno do ângulo de azimute de 90°, que tem menos energia em ângulos de reflexão afastados porque o deslocamento de linha transversa máxima (azimute de superfície de 90°) é somente metade do deslocamento em linha máximo (azimute de superfície de 0o).

De acordo com uma concretização exemplíficativa ilustrada na Figura 23, um método para calcular as CIG inclui uma etapa 2300 de calcu- lar uma fonte de campo de ondas (baseada na equação (20)) e uma etapa 2302 de calcular um campo de ondas de receptor (baseada na equação (21)). Os campos de ondas gerados nas etapas 2300 e 2302 são CIGs de domínio de ângulo de amplitude real de 3D provenientes da RTM. Esta for- mulação é numericamente estável, adapta-se a geometrias de aquisição di- ferentes, proporciona entrada de alta qualidade para tomografia, e propor- ciona amplitude de migração favorável para a análise de amplitude VS ângu- lo.

Para reduzir o custo de computação, janelas locais são aplica- das na etapa 2304. Não obstante, esta etapa pode introduzir efeitos de limi- te indesejáveis. Desta forma, na etapa 2306 é aplicada uma sobre amostra- gem no domínio de número de ondas. Muito embora reduza os efeitos de limite, esta etapa aumenta o custo computacional. Um algoritmo ALFT é a- plicado na etapa 2308 a cada um dos resultados das etapas 2300 e 2302. A ALFT sem etapas adicionais pode aumentar o custo computacional. Desta forma, na etapa 2310 somente são selecionadas aquelas contribuições que são dotadas da energia mais alta, o que reduz o custo computacional tal como discutido anteriormente. Na etapa 2310 a Fonte e campos de ondas de receptor no domínio k são are obtidos e estes campos de ondas têm uma representação esparsa devido à ALFT. Desta forma, na etapa 2312, uma condição de formação de imagem do tipo correlação (por exemplo, equação (24)) é aplicada ao domínio de número de ondas e os ângulos Θ e cp são cal- culados, por exemplo, com base na equação (1). As CIG são determinadas na etapa 2314. Não obstante, uma vez que as CIG estão no domínio k, uma FFT inversa é aplicada na etapa 2316 para determinar as CIG no domínio de espaço tal como ilustrado na etapa 2318. O método discutido anteriormente pode ser implementado em uma concretização exemplificativa em um aparelho de processamento 2400 tal como ilustrado na Figura 24. O dispositivo de processamento 2400 pode ser configurado especificamente para calcular convergências de imagens comuns de domínio de ângulo (ADCIGs) para o uso na geração de uma i- magem final de uma subsuperfície da terra. O dispositivo de processamento 2400 pode ser, por exemplo, um computador, um processador, um servidor, ou uma rede de computadores, processadores ou servidores. O dispositivo de processamento 2400 pode incluir uma interface 2402 configurada para receber dados vindos do exterior, por exemplo, dados sísmicos. A interface 2402 pode acomodar conexão da Internet, conexão por satélite, teclados e assemelhados. O dispositivo de processamento 2400 também inclui um processador 2404 conectado à interface 2402 e configurado para executar uma ou mais ou todas as etapas discutidas com relação à 23. O dispositivo de processamento 2400 pode ser dotado de circuitos dedicados para cada etapa da Figura 23 ou o processador 2404 pode ser configurado com softwa- re para realizar todas as etapas ilustradas na Figura 23. O dispositivo de processamento 2400 também pode incluir um display 2406 para exibir uma imagem final calculada por meio do processador 2404. Vários dos dados usados para calcular a imagem final podem ser armazenados em um dispo- sitivo de armazenamento 2408 o qual é conectado ao processador 2404.

Também poderão estar presentes outros componentes conhecidos de um computador, servidor ou processador.

As concretizações exemplificativas expostas proporcionam um sistema e um método para processar com mais precisão e rapidez dados relacionados com uma subsuperfície. Deverá ser compreendido que esta descrição não tem por finalidade limitar a invenção. Pelo contrário, as con- cretizações exemplificativas destinam-se a abranger alternativas, modifica- ções e equivalentes, que estão incluídos no espírito e escopo da invenção tais como definidos pelas reivindicações em anexo. Além disso, na descri- ção detalhada das concretizações exemplificativas, encontram-se expostos numerosos detalhes específicos com a finalidade de proporcionar uma com- preensão ampla da invenção reivindicada. Não obstante, aquele versado na técnica irá compreender que várias concretizações podem ser praticadas sem esses detalhes específicos.

Muito embora os aspectos e elementos das presentes concreti- zações exemplificativas estejam descritos nas concretizações em combina- ções particulares, cada aspecto ou elemento pode ser usado isoladamente sem os outros aspectos e elementos das concretizações ou nas várias com- binações com ou sem outros aspectos e elementos expostos neste contexto.

Esta descrição por escrito utiliza exemplos da matéria em apreço expostos para permitir a qualquer pessoa versada na técnica pôr os mesmos em prática, incluindo a execução e utilização de quaisquer dispositivos ou sistemas e a realização de quaisquer métodos incorporados. O escopo pa- tenteável da matéria em apreço é definido por meio das reivindicações, e pode incluir outros exemplos que ocorrerão àqueles versados na técnica.

Esses outros exemplos destinam-se a ficar incluídos dentro do escopo das reivindicações.

Claims (10)

1. Método para calcular, em um dispositivo de computação, con- vergências de imagens comuns de domínio de ângulos (ADCIG) para o uso na geração de uma imagem final de uma subsuperfície da terra, o método compreendendo: calcular (2300) com propagação para a frente de campos de on- das, em um domínio de tempo espaço, uma fonte de campo de ondas Pf de uma fonte sísmica; calcular (2302) com propagação para trás de campos de ondas, no domínio de espaço de tempo, um campo de ondas de receptor pB de um receptor sísmico configurado para detectar uma onda sísmica gerada pela fonte sísmica, em que a onda sísmica emitida pela fonte sísmica entra na subsuperfície da terra e é refletida em um refletor da subsuperfície antes de chegar ao receptor sísmico; aplicar (2308) um algoritmo de transformação de Fourier anti- vazamento (ALFT) para transformar o campo de ondas de fonte Pf para um domínio de número de ondas; aplicar (2308) o algoritmo de ALFT ao campo de ondas de re- ceptor para transformar o campo de ondas de receptor no domínio de núme- ro de ondas; determinar (2312) uma condição de formação de imagem para a fonte de ALFT e campos de ondas de receptor no domínio de número de ondas; computar (2312) um ângulo de reflexão Θ e um ângulo de azimu- te φ do campo de ondas de fonte pf e campo de ondas de receptor pB no domínio de número de ondas; calcular (2314) as ADCIG no domínio de número de ondas; e aplicar (2316) uma transformação de Fourier rápida inversa (FFT) para determinar as ADCIG no domínio de espaço.

2. Método, de acordo com a reivindicação 1, que compreende ainda: calcular o campo de ondas de fonte Pf como uma propagação para a frente de uma assinatura de fonte.

3. Método, de acordo com a reivindicação 2, em que o campo de ondas de fonte pF é calculado com base nas equações: em que v é a velocidade de pF, V2 é o operador de Laplacian, x é uma posição onde o operador de Laplacian é aplicado ao campo de ondas de fonte pF, t é o tempo em que o campo de ondas de fonte pF está presente na posição x, z é um componente de profundidade da posição x, ^(x-s)/(/) é a assinatura de fonte, δ é uma função de Dirac, e s é uma localização da fonte sísmica.

4. Método, de acordo com a reivindicação 1, que compreende ainda: calcular o campo de ondas de receptor pB como uma propaga- ção para trás dos dados registrados pelo receptor sísmico.

5. Método, de acordo com a reivindicação 4, em que o campo de ondas de receptor pB é calculado com base nas eauacões: em que v é a velocidade de pB, V2 é o operador de Laplacian, x é uma posição onde o operador de Laplacian é aplicado ao campo de ondas Pr de refletor, t é o tempo no qual o campo de ondas de receptor pB está presente na posição x, z é uma profundidade da posição x, dGobs(x,y,$‘,t) é constituído pelos dados registrados, e s é a localização da fonte sísmica.

6. Método, de acordo com a reivindicação 1, que compreende ainda: determinar uma equação assintótica de migração de amplitude real tanto no ângulo de reflexão Θ quanto no ângulo de azimute cp, em que a equação é: em que R é a refletividade da subsuperfície, x é a posição onde a refletividade é calculada, s é uma posição da fonte sísmica, r é uma posi- ção do receptor sísmico, SGohs ê compreendido pelos dados registrados pelo receptor sísmico, δ é uma função de Dirac, ω é uma freqüência dos campos de ondas, A é um produto entre (i) uma amplitude de raio a partir da locali- zação s para a localização x e (ii) uma amplitude de raio desde a localização x até à localização r, p] é relacionado com pF e p\ é relacionado com pB.

7. Método, de acordo com a reivindicação 6, que compreende ainda: transformar a equação assintótica de migração de amplitude real em uma equação de migração de tempo inverso de amplitude real (RTM) que é dada por meio de em que as convergências de imagem comuns (CIG) são dadas por Έ(χ,θ0,φ0) = R(x, 0o,<po)v(x)/2cos0o, pF =G*(x,s)pl , e com G sendo uma função de Green.

8. Método, de acordo com a reivindicação 7, em que as ADCIG são dadas por: com k sendo o número de ondas.

9. Método, de acordo com a reivindicação 8, em que as ADCIG são reescritas como mm haco na mnrtinãn Ηω r~"TiaÇãO (jg imagem expressa no domínio de número de ondas.

10. Dispositivo de processamento (2400) configurado para calcu- lar convergências de imagens comuns de domínio de ângulo (ADCIG) para o uso na geração de uma imagem final de uma subsuperfície da terra, com o dispositivo de processamento que compreende: uma interface (2402) configurada para receber dados sísmicos; e um processador (2404) conectado à interface (2402) e configu- rado para, calcular (2300) com propagação para a frente de campos de on- das, em um domínio de espaço de tempo, uma fonte de campo de ondas pF de uma fonte sísmica; calcular (2302) com propagação para trás de campos de ondas, no domínio de espaço de tempo, um campo de ondas de receptor pB de um receptor sísmico configurado para detectar uma onda sísmica gerada pela fonte sísmica, em que a onda sísmica emitida pela fonte sísmica entra na subsuperfície da terra e é refletida em um refletor da subsuperfície antes de chegar ao receptor sísmico; aplicar (2308) um algoritmo de transformação de Fourier anti- vazamento (ALFT) para transformar o campo de ondas de fonte pF para um domínio de número de ondas; aplicar (2308) o algoritmo de ALFT ao campo de ondas de re- ceptor para transformar o campo de ondas de receptor no domínio de núme- ro de ondas; determinar (2312) uma condição de formação de imagem para a fonte de ALFT e campos de ondas de receptor no domínio de número de ondas; computar (2312) um ângulo de reflexão Θ e um ângulo de azimu- te φ do campo de ondas de fonte pF e campo de ondas de receptor pB no domínio de número de ondas; calcular (2314) as ADCIG no domínio de número de ondas; e aplicar (2316) uma transformação de Fourier rápida inversa (FFT) para determinar as ADCIG no domínio de espaço.