BR112016016088B1

BR112016016088B1 - Método implementado por computador de migração de dados sísmicos

Info

Publication number: BR112016016088B1
Application number: BR112016016088-6A
Authority: BR
Inventors: Sheng Xu; Hongbo Zhou
Original assignee: Statoil Petroleum As
Priority date: 2014-01-10
Filing date: 2015-01-09
Publication date: 2023-02-23
Also published as: EA201691406A1; US20200271802A1; US11092707B2; US10705233B2; CN106133556B; CA2936326A1; AU2015205510A1; EA031826B1; MX2016009033A; BR112016016088A2; AU2015205510B2; US20160334527A1; WO2015104059A1; CN106133556A; WO2015104386A1

Abstract

DETERMINAÇÃO DE UMA COMPONENTE DE UM CAMPO DE ONDA. A presente invenção refere-se a modalidades descritas com relação a um método para determinar um campo de onda em uma subsuperfície anisotrópica da Terra. O método inclui resolver numericamente uma equação de onda de modo de onda única quase acústica com base em parâmetros anisotrópicos espacialmente variados, para determinar o campo de onda na subsuperfície anisotrópica.

Description

CAMPO DA INVENÇÃO

[0001] A presente invenção refere-se à determinação de um campo de onda ou uma componente de um campo de onda na subsuperfície da terra, em particular uma subsuperfície anisotrópica.

FUNDAMENTOS DA INVENÇÃO

[0002] A propagação de uma onda através de um meio pode ser descrita por uma “equação de onda”. A equação de onda pode, tipicamente, ser resolvida para determinar a magnitude ou outra propriedade da onda como uma função do tempo e do espaço.

[0003] Levantamentos geofísicos podem ser realizados para explorar a subsuperfície da Terra. Por exemplo, levantamentos sísmicos ou eletromagnéticos podem ser realizados se uma onda acústica ou eletromagnética (EM) é transmitida a partir de uma fonte para a subsuperfície e detectores são utilizados para detectar uma resposta. Com base na resposta detectada, as propriedades da terra podem ser determinadas. Por exemplo, a resposta pode compreender dados de séries temporais com características associadas com interfaces na subsuperfície.

[0004] No levantamento de reflexão sísmica, eventos de grande amplitude em dados sísmicos podem estar relacionados a interfaces através das quais as propriedades das rochas mudam abruptamente e refletem uma quantidade significativa de energia de volta para os detectores. Por conseguinte, as altas amplitudes nos dados podem ser indicativas de interfaces de subsuperfície.

[0005] No processamento típico de dados sísmicos, o tempo entre a transmissão da onda sísmica e a detecção do evento de alta amplitude é frequentemente considerado como sendo o tempo de viagem de uma onda sísmica para a interface e de volta, devido a uma reflexão no ponto médio entre a fonte e o receptor. O tempo de viagem é, portanto, considerado tipicamente um aviso para a profundidade, e a plotagem de dados sísmicos como uma função do tempo de viagem pode, portanto, indicar a relação de profundidade de características de amplitude e, portanto, de interfaces de subsuperfície.

[0006] Na realidade, no entanto, a velocidade sísmica das camadas da subsuperfície através das quais a onda passou não é constante (tipicamente aumenta com a profundidade), e, por conseguinte, a fim de se obter uma impressão mais realista da localização das interfaces, os dados são convertidos em profundidade, isto é, representados graficamente como uma função da profundidade, com base em um perfil de velocidade ou modelo de velocidade predeterminado.

[0007] Os dados para os pares de fonte-receptor são representados graficamente em um ponto médio entre eles. Os dados de todos esses pares podem, portanto, aparecer em um registro de tempo ou profundidade nesse ponto médio, com eventos de amplitude que aparecem como se derivassem de um arranjo fonte-receptor de incidência normal. Enquanto isso pode fornecer uma representação adequada da estrutura da Terra onde as camadas são horizontais e planas, está claro que uma onda sísmica pode ser refratada ao longo do seu caminho de propagação, ser refletida a partir de interfaces inclinadas e interagir com diferentes camadas nas direções incidente e refletida. Por conseguinte, os registros de dados convertidos de profundidade e dados de tempo de viagem podem sugerir uma falsa geometria. Por exemplo, um evento de alta amplitude pode ocorrer nos dados a uma posição de profundidade particular onde, na realidade, não existe nenhum refletor nessa posição de profundidade.

[0008] Este é um problema conhecido no processamento de dados sísmicos e na interpretação dos dados de reflexão sísmica, particularmente em estruturas geologicamente complexas. A correção dos dados pode ser tentada por um processo de migração. A migração sísmica serve para reposicionar os dados sísmicos para suas posições representativas, geometricamente fiéis, no conjunto de dados. Algoritmos de migração foram desenvolvidos para a migração de dados sísmicos. É útil compreender a propagação de onda da onda sísmica através da subsuperfície, a fim de realizar a migração. Em termos simples, traçando a rota que a onda tomou da fonte para o refletor e de volta para a superfície, é possível identificar os receptores nos quais a energia associada com esse refletor chegou e o tempo de chegada dessa energia. A fim de fazer isso, é necessário um modelo de velocidade da subsuperfície. Este pode ser obtido a partir de um pacote de construção de modelo e pode ser baseado em dados de poço ou obtidos de outras maneiras, por exemplo, estimados, talvez baseadas no conhecimento geológico regional.

[0009] Outro fator complicador é a anisotropia da subsuperfície. Por exemplo, determinados tipos de rochas podem ter uma forte direcionalidade de propagação. As velocidades sísmicas em diferentes direções (por exemplo, ortogonais) podem ser significativamente diferentes, e podem diferir entre as camadas.

[0010] Técnicas de migração de campo de onda total foram desenvolvidas para modelar o desenvolvimento da onda no tempo e em posições espaçadas espacialmente em 3D através de uma região da subsuperfície. Tais técnicas podem envolver resolver uma equação de onda do tipo mencionado acima, para obter o campo de onda em diferentes localizações, em relação ao tempo em relação a uma emissão de sinal de fonte em diferentes localizações.

[0011] Esses algoritmos de migração têm como uma entrada dados sísmicos padrão. Os dados sísmicos poderiam ser organizados como ganhos de disparo comuns, ou em algum outro domínio, tal como onda plana ou domínios de disparo atrasados. Os dados sísmicos podem ser fornecidos em um domínio transformado (em comparação com os dados sísmicos originalmente adquiridos), por exemplo, um domínio combinado de fase linear para localização de disparo para a migração de disparo retardado, em que a combinação dos dados é realizada por transformada tau- p para a localização do disparo. A transformada tau-p pode também ser realizada para a localização do receptor para realizar uma migração de onda plana. A solução de campo de onda da equação de onda fornece um conjunto de valores complexos em cada localização/ponto alvo, fornecendo, por exemplo, amplitude, fase, número de onda, etc. Esses valores são as componentes de entrada cruciais para os algoritmos de migração.

[0012] Muitas técnicas de equações de onda diferentes têm sido desenvolvidas, e é de interesse desenvolver uma técnica de equação de onda para utilização em migração sísmica no caso de uma subsuperfície anisotrópica. Algumas técnicas conhecidas de migração são agora descritas em mais detalhes.

[0013] Hoje em dia, modelos anisotrópicos sísmicos precisos são reconhecidos como pontos-chave na captação de imagens de áreas muito desafiadoras para estruturas complexas, porque descrevem a discrepância da velocidade de propagação direcional das ondas sísmicas, que é muito mais realista para a estrutura geológica sobrecarregada tal como sedimentos de leito com fraturas. A isotropia transversal inclinada (TTI) ou anisotropia ortorrômbica inclinada (TOA) são geralmente necessárias para a captação de imagem de tais estruturas anisotrópicas complicadas (Zhang e Zhang, 2008; Zhang e Zhang, 2011; Fowler e King, 2011). Dado um modelo anisotrópica preciso, algoritmos de migração são destinados a fornecer as imagens de subsuperfície livre de artefatos. Isto pode exigir aplicações de algoritmos de migração de mínimos quadrados (Lambaré e outros, 1992; Nemeth e outros, 1999) e pode tornar-se ainda mais complicado quando se trabalha em captação de imagem do modelo elástico (Jin e outros, 1992). Os algoritmos de migração baseados em integrada atuais, tais como migrações de Kirchhoff, são conhecidos por serem muito eficazes à medida que o traçado de raios pode ser facilmente adaptado para a equação Eikonal anisotrópica (Han e Xu, 2012). Estes algoritmos são desenvolvidos para um único modo de ondas, tais como a onda quase-P, que permite que os algoritmos de migração sejam linearizados, de modo a reduzir o custo computacional. No entanto, para os algoritmos baseados em equação de onda, uma simulação eficiente e precisa de uma onda de único modo, tal como onda quase-P ou quase-S, é essencial para o desenvolvimento de algoritmos de migração de sucesso, tais como a migração reversa no tempo (RTM) (Baysal e outros, 1983; McMechan, 1983; Whitmore, 1983).

[0014] Historicamente, existem maneiras diferentes para simulação numérica de modo único de ondas. A primeira delas consiste em resolver as equações de onda elástica completas e, em seguida, dividir os campos de onda para separar a onda quase-P para posterior análise. A separação de campo de onda pode ser eficaz para o caso isotrópico (Sun e outros, 2004), mas não é uma tarefa fácil para os casos anisotrópicos heterogêneos e pode ser extremamente computacionalmente ineficiente (Yan e Sava, 2009; Cheng e Fomel, 2013).

[0015] Outra maneira de reduzir o custo computacional enquanto mantendo a propagação da onda anisotrópica isotrópica transversal (TI) é aplicar uma aproximação acústica (Alkhalifah, 2000) à equação de TI, definindo a velocidade de onda de cisalhamento ao longo do eixo de simetria para zero. Isto leva a uma equação diferencial de quarta ordem escalar. No entanto, Alkhalifah (2000) não propõe qualquer técnica numérica viável para obter o campo de onda a partir desta equação diferencial de quarta ordem. Zhou e outros (2006) decompõem esta equação diferencial de quarta ordem em um sistema acoplado de equações diferenciais de segunda ordem 2x2. Isso resulta em um esquema mais computacionalmente eficiente do que as equações elásticas originais. No entanto, as aplicações para o meio TTI com o eixo simétrico variável TTI, especialmente com a existência de mudanças abruptas no eixo simétrico TTI, demonstram que este sistema não é numericamente estável: surgem as instabilidades fracas e os ruídos crescem linearmente no tempo (Liu e outros, 2009, Zhang e outros, 2011). Estas instabilidades foram bem analisadas e soluções foram propostas (Bakker e Duveneck, 2011; Zhang e outros, 2011, Bube e outros, 2012.). Embora as instabilidades pudessem ser resolvidas, definir zero para a velocidade de onda de cisalhamento ao longo do eixo de simetria não exclui a existência das ondas de pseudo-cisalhamento, que são as soluções intrínsecas deste sistema acoplado. As pessoas têm tentado de diversas maneiras atenuá-las (Zhang e outros, 2009; Guan e outros, 2011); ou simplesmente ignorar a presença desses artefatos de onda de cisalhamento como práticas atuais da indústria de TTI RTM, com a expectativa de que fossem anulados por condições de imagem e pilha de migração.

[0016] Outra maneira é computar diretamente a equação da onda P quase- acústica desacoplada. Dados os parâmetros anisotrópicos constantes ε e δ, a equação desacoplada poderia ser resolvida com o algoritmo pseudo espectral (Etgen e Brandsberg-Dahl, 2009). Portanto, permite-nos tomar um conjunto de soluções de diferentes constantes ε e δ, prosseguindo com um esquema de interpolação, para resolver numericamente a onda quase-P com os parâmetros anisotrópicos que estão mudando levemente no domínio espacial (Chu e outros, 2013). Para meios VTI (isotropia transversal com um eixo de simetria vertical), o esquema de interpolação é bidimensional com relação a ε e δ. Mas para meios TTI, se o eixo de simetria também varia espacialmente, duas dimensões adicionais são necessárias para a interpolação. Esse método de interpolação mais pseudo-espectral pode precisar computar um conjunto enorme de combinações de parâmetros anisotrópicos e pode ser muito ineficaz quando é aplicado a um modelo muito complicado. Quando as complexidades do modelo anisotrópico são leves, a equação desacoplada pode ser resolvida pelo método separável aproximadamente otimizado (Liu e outros, 2009), o que pode ser considerado como um caso especial de aproximação de baixa classificação (Fomel e outros, 2012). O custo computacional pode ser menor do que método de interpolação mais pseudo-espectral, mas ainda muito ineficiente e impreciso.

SUMÁRIO DA INVENÇÃO

[0017] Vários aspectos da invenção são fornecidos como apresentados nas reivindicações em anexo.

[0018] Qualquer aspecto pode ter outras características, tal como descrito em relação a qualquer outro aspecto e as características entre aspectos podem ser combinadas e intercambiadas.

[0019] As modalidades da invenção fornecem, vantajosamente, as técnicas para calcular numericamente o campo de onda de uma onda se propagando em um meio anisotrópico. Em particular, as técnicas são vantajosas em termos de custo numérico ou computacional reduzido e simplicidade de implementação.

[0020] A tecnologia existente geralmente usa sistemas de equação acoplada que descrevem a propagação de onda e fornecem soluções de campo de onda para vários modos de onda juntos. Tipicamente, sistemas acoplados para os modos de onda quase-P e quase-S são utilizados. Tal sistema tem, intrinsecamente, dois valores próprios associados com os modos de onda quase-P e quase-S, respectivamente. Em técnicas de equação acoplada existentes, artefatos e ruídos nas soluções de campo de onda aparecem devido à inclusão do modo de onda quase-S. As soluções de equações acopladas geralmente têm alto custo numérico.

[0021] Têm faltado técnicas para calcular numericamente o campo de onda de uma onda sísmica se propagando em regiões anisotrópicas da Terra.

BREVE DESCRIÇÃO DOS DESENHOS

[0022] Serão descritas agora, a título de exemplo apenas, as modalidades da invenção com relação aos desenhos em anexo, nos quais:

[0023] A Figura 1 é uma representação em fluxograma de um método para determinar um campo de onda de acordo com uma modalidade da invenção.

[0024] A Figura 2 é uma representação esquemática de um dispositivo de computador para utilização na determinação do campo de onda de acordo com uma modalidade da invenção.

[0025] A Figura 3 é uma representação esquemática do aparelho para adquirir dados.

[0026] A Figura 4 é um gráfico de contorno de magnitudes de campo de onda para uma solução de campo de onda de modelo VTI.

[0027] A Figura 5 é um gráfico de contorno de magnitudes de um operador utilizado para calcular o campo de onda do modelo da Figura 4.

[0028] A Figura 6 é um gráfico de contorno de magnitudes de um operador composto S| V U no campo de onda das Figuras 4 e 5.

[0029] A Figura 7 é um gráfico de contorno de magnitudes de campo de onda para uma solução de campo de onda de modelo TTI.

[0030] A Figura 8 é um gráfico de contorno de magnitudes de um operador utilizado para calcular o campo de onda do modelo da Figura 7.

[0031] A Figura 9 é um gráfico de contorno de magnitudes de um operador composto S| V U no campo de onda das Figuras 7 e 8.

[0032] A Figura 10 é um gráfico de contorno do modelo de velocidade para um modelo de sal SEAM TTI.

[0033] A Figura 11 é um gráfico de contorno de magnitudes de campo de onda para o modelo da Figura 10 no tempo t = 1,6s.

[0034] A Figura 12 é um gráfico de contorno de magnitudes de campo de onda para o modelo da Figura 10 no tempo t = 2,4s.

[0035] A Figura 13 ilustra três curvas de dispersão.

[0036] A Figura 14 ilustra a amplitude dos termos de correção assintótica.

[0037] A Figura 15 ilustra dois campos de onda.

[0038] A Figura 16 ilustra um cálculo executado com um modelo elíptico.

DESCRIÇÃO DETALHADA DA INVENÇÃO

[0039] Modelo de equação de onda

[0040] Primeiro propõe-se uma equação de onda quase-P pura cinematicamente precisa para modelos anisotrópicos complexos. Esta equação é obtida pela decomposição do operador pseudo-diferencial original da Alkhalifah (2000) em dois termos numéricos solucionáveis: um operador diferencial e um operador escalar. A equação de onda quase-P nova resultante separou os parâmetros do modelo anisotrópico do operador pseudo-diferencial espacial e, portanto, pode ser resolvida com esquemas numéricos tradicionais sem sacrificar a sua precisão, ou seja, nenhuma aproximação foi aplicada e, portanto, mantém a mesma relação de dispersão da equação de onda pseudo-diferencial original. Os operadores diferenciais na equação são projetados para serem auto-adjuntos e, portanto, podem conservar a energia quando a onda se propaga. Em seguida, discute-se a implementação numérica do algoritmo proposto. A equação de onda quase-P proposta tem uma forma similar a uma equação de onda acústica e pode ser implementada de forma similar. Nesta implementação atual, usa-se o algoritmo pseudo-espectral para as derivadas espaciais e o esquema de diferenças finitas de precisão de segunda ordem para a derivada temporal. Extensões do novo algoritmo para meios TTI e TOA são diretas e foram brevemente mencionadas nos textos.

[0041] Valida-se esta abordagem com exemplos das respostas de impulso dos modelos VTI/TTI e ilustra-se as funcionalidades dos termos na nova equação de onda quase-P pura. Por fim, demonstra-se a eficácia e robustez desta equação quase-P com um exemplo das propagações de onda quase-P em um modelo de sal SEAM TTI.

[0042] Decomposição de operador pseudo-diferencial

[0043] Em seguida, para simplificar trabalha-se no primeiro caso VTI. Isto faz com que se resolva a seguinte equação pseudo-diferencial escalar (Alkhalifah, 2000):

ao longo do eixo de simetria; ε e δ são os parâmetros anisotrópicos de Thomsen (1986), como definido por Thomsen (1986). As definições destes parâmetros, como fornecidas por Thomsen (1986, ver abaixo a seção de referência) são incorporadas na presente descrição por referência. Nota-se que todos estes parâmetros são espacialmente variados. A Eq. 1 regula a propagação da onda quase-P, embora a amplitude de sua solução possa ser bastante diferente de uma onda quase-P real elástica. Isto ocorre devido ao desacoplamento das equações elásticas completas após o que toda as interações de diferentes modos de onda são omitidas, incluindo a conversão para e a partir de outros modos de onda. No entanto, A Eq. 1 compartilha a mesma dispersão com essa onda quase- P elástica, indicando que a sua fase é precisa em comparação com a onda quase-P elástica. Embora a Eq. 1 seja muito promissora para a propagação de onda de modo único, infelizmente, ela representa uma equação de operador pseudo-diferencial que não pode ser resolvida com esquemas numéricos de diferença finita tradicionais. Alkhalifah (2000) não resolve numericamente a Equação 1 em meios complexos.

[0044] A formulação desse operador pseudo-diferencial sugere uma solução numérica na qual as componentes de gradiente espacial do campo de onda atual são determinadas primeiro, formando a combinação dos componentes de gradiente dentro da raiz quadrada, e, em seguida, tomando a raiz quadrada dessa combinação. Isso não funciona em esquemas tradicionais. Reescreve-se a equação, e substitui-se o operador pseudo-diferencial por dois operadores que não são mais operadores pseudo-diferenciais.

[0045] De modo a descrever esta abordagem, inicia-se com a sua relação de dispersão correspondente da Eq. 1, que é representada como segue:

Aqui a é a frequência angular. O vetor de número de onda espacial k é habitualmente definido como kR = (kx,ky,kz ). Portanto, kh é o número de onda horizontal com kh2 =kx2 +ky2 , kz é o número de onda vertical. A Eq. 2 também pode ser reescrita como:

Onde nR= (nx,ny,nz) é o vetor unitário de direção de fase e é definido como

[0046] Define-se um operador escalar auxiliar S como

[0047] Agora, a Eq. 3 torna-se

[0048] No domínio espacial, o operador -k2 é expresso como V*V onde o símbolo V^ denota divergência, e V denota gradiente. Transformar a Eq. 6 do domínio da frequência-número de onda de volta para o domínio do tempo-espaço resulta na seguinte equação: d2,,

[0049] Esta equação representa as equações diferenciais parciais correspondentes da Eq. 6 no domínio de tempo-espaço. Isto define um modelo para a propagação da onda quase-P em um meio anisotrópico, tal como na subsuperfície da Terra, solucionável e fácil de implementar e métodos numéricos eficientes para fornecer um campo de onda no espaço e no tempo. Não é mais uma equação pseudo- diferencial porque separou-se os termos anisotrópicos das derivadas. Mas esta equação é uma equação de onda não linear uma vez que o operador escalar S depende também da solução dos campos de onda. Este operador escalar pode ser calculado no domínio do espaço à medida que o vetor nR tem o seu significado físico: direção de fase da frente de onda. Como a Eq. 5 não é sensível ao sinal das componentes de nR (uma vez que todas as componentes têm número par de energia), u . Vu . . R pode-se apenas utilizar as componentes de ,—f para substituir as componentes de n ' ' ' Vu| ' ' na Eq. 5. O operador escalar S é crucial nesta abordagem; ele realmente desempenha o papel da relação de dispersão por toda a propagação da onda. Ele controla a velocidade de propagação da onda quase-P (um operador de propagação anisotrópica), ele depende não somente da direção de fase da propagação em cada ponto espacial, mas também dos parâmetros anisotrópicos em cada localização espacial individual. E a computação deste operador precisa apenas do parâmetro anisotrópico que varia espacialmente, e o gradiente do campo de onda atual.

[0050] Pode notar-se que a propagação de onda é calculada em intervalos de tempo. A fim de computar o próximo intervalo de tempo do campo de onda, é necessário o campo de onda atual. Um campo de onda de partida é definido, fornecendo uma condição limite ou ponto de partida para os cálculos. Uma ondeleta de origem predefinida emitida a uma localização de origem pode ser definida e utilizada para permitir que este campo de onda de partida seja determinado.

[0051] O modelo de velocidade é predeterminado, por exemplo, estimado por outros algoritmos e/ou fornecido em um pacote de construção de modelo. O pacote de construção de modelo pode incluir muitos algoritmos tais como inversão de tomografia, inversão de onda completa, interpretação de sal, etc. A velocidade é uma entrada necessária e deve ser conhecida por algoritmos de migração. Os parâmetros de anisotropia são partes do modelo para a propagação de onda (entrada), e também contidos no modelo de velocidade. Eles poderiam ser estimados pelo pacote de construção de modelo.

[0052] O operador diferencial na Eq. 7 é um operador auto-adjunto e, portanto, conserva a energia. Isso deve garantir a estabilidade das propagações de onda mesmo para os casos com alterações espaciais abruptas dos parâmetros de modelo e anisotrópicos.

[0053] A generalização de VTI para TTI ou TOA é muito simples nesta abordagem. Para TTI precisa-se apenas projetar o vetor de gradiente de campo de onda Vu para as coordenadas locais, em que o eixo z local é o eixo de simetria da anisotropia, e depois aplicar exatamente o mesmo procedimento para a Eq. 5. Uma análise similar pode ser generalizada aos meios TOA.

[0054] Uma solução alternativa de Eq. 6 pode ser implementada com a equação descrita como segue:

[0055] Aqui, o operador pseudo-diferencial na Eq. 1 é decomposto em dois operadores: um operador escalar e um operador de Laplace. Esta equação utiliza diretamente o operador de Laplace e tem uma forma similar a uma equação de onda acústica.

[0056] Decomposição em termos esféricos

[0057] O cálculo numérico do operador escalar S no domínio de número de onda exige que os parâmetros anisotrópicos na equação 7 sejam constantes, e pode, portanto, ser calculado no domínio do espaço com a direção de fase aproximada assintoticamente: n = V u /| V u\

[0058] Para analisar teoricamente o efeito do termo de aproximação assintótica, reescreve-se o operador escalar na equação 7 como:

[0059] Usando esta expressão para o operador escalar, a equação de onda se torna

[0060] O primeiro termo do lado direito na Equação 10 é a equação de onda de fundo. Pode ser considerada como um operador diferencial - o operador de Laplace, que não contém uma aproximação. O segundo termo no lado direito da Equação 10 pode ser considerado como um termo de correção. O cálculo de ΔS depende das direções de propagação de onda, o que é uma aproximação assintótica. A Equação 10 pode ser considerada como uma decomposição esférica da equação de onda original 1.

[0061] Na Figura 13, a curva de dispersão 131 corresponde a uma solução da Equação 1, a curva de dispersão 132 é uma curva de dispersão de fundo correspondente ao primeiro termo do lado direito da Equação 10. A curva de dispersão 133 será discutida em relação a um modelo elíptico abaixo. A dispersão de fundo 132 é ilustrada para os parâmetros anisotrópicos constantes ε = 0,25 e δ = 0,1. A correção assintótica é destinada a corrigir a fase a partir da curva 132 até a curva 131.

[0062] Na Figura 14, a curva 141 mostra a magnitude da correção assintótica ΔS entre os ângulos 0 a 360°. Nota-se que a correção é relativamente grande em comparação com a dispersão de fundo: para este exemplo, a correção máxima é aproximadamente 50% da dispersão de fundo. Por conseguinte, esta abordagem exige uma estimativa muito precisa da direção do campo de onda. Uma amostragem espacial lateral de computações maiores do que a amostragem de Nyquist pode ser utilizada. Portanto, a parte de alto número de onda dos campos de onda é serrilhada e pode introduzir erros no cálculo do vetor de direção na faixa de número de onda inferior de campos de onda, o que torna esta abordagem vulnerável aos erros de direção.

[0063] Decomposição em termos elípticos

[0064] A fim de melhorar a tolerância dos erros numéricos do vetor direcional no algoritmo, propõe-se, em vez de uma decomposição esférica, uma decomposição elíptica como definida na Equação 11 abaixo. Nesta equação, decompõe-se ainda o operador pseudo-diferencial em dois operadores: um operador diferencial e um operador escalar. No entanto, o operador de Laplace na decomposição original é substituído por um operador diferencial elíptico, enquanto o operador escalar é também correspondentemente modificado para assegurar a fase precisa da propagação de onda. O objetivo da nova decomposição é reduzir a magnitude do termo assintótico. Esta decomposição é:

onde Se é o escalar elíptico. O termo

pode ser escrito como

e interpretado como os tempos de operador diferencial u e este operador diferencial (o termo dentro dos grandes parêntesis) é um operador diferencial elíptico. Para analisar melhor o termo assintótico na Equação 5, reescreve-se a mesma explicitamente como (Equação 12): d

[0065] A Figura 13 mostra as comparações das curvas de dispersão entre a abordagem esférica e a abordagem elíptica. Para os parâmetros anisotrópicos constantes ε = 0,25 e δ = 0,1, a curva 133 mostra a curva de dispersão do primeiro termo do lado direito da Equação 12. Comparando esta curva com a abordagem de decomposição esférica, nota-se que o fundo de decomposição elíptica está muito mais próximo da solução exata desejada ilustrada como a curva 131. Na Figura 14, a linha 142 mostra a magnitude da correção assintótica ΔSe em 0 ~ 360°. A magnitude máxima de ΔSe é 0,068, que é 7 vezes menor do que a decomposição esférica. Portanto, a decomposição elíptica tem tolerância muito melhor aos erros de direção.

[0066] Primeiro, demonstra-se os efeitos deste algoritmo com um exemplo de um simples impulso na Figura 15. Este exemplo é um modelo TTI simples (isotrópico transversalmente inclinado), que é homogêneo com a velocidade vertical definida como 2000 m/s e os parâmetros anisotrópicos de Thomsen (1986) como ε = 0,24 e δ = 0,1 com um ângulo inclinado de 30° e um ângulo de azimute de 135°. A ondeleta de origem para este exemplo é uma ondeleta de Ricker com uma frequência máxima de 24 Hz. A grade de computação é um cubo 3D com comprimentos de 6,0 km e a amostragem espacial como 15 m em todas as 3 direções. A fonte pontual é colocada no meio da grade. A Figura 15a) traça uma fatia sísmica 2D, localizada no meio da direção Y, do instantâneo de campo de onda 3D no tempo t = 1,0 s a partir da abordagem de decomposição esférica; e a Fig 15b) traça a mesma com abordagem de decomposição elíptica. Ambas as abordagens geram apenas onda pseudo-P e não existe nenhuma onda de cisalhamento. Ambos os campos de onda produzem a mesma fase de propagação, mas a abordagem de decomposição elíptica fornece uma amplitude mais equilibrada. Observa-se ainda que os custos numéricos para decomposição esférica e para a decomposição elíptica são quase idênticos.

[0067] O segundo exemplo é um teste de migração com o modelo SEG SEAM. Selecionou-se uma linha de disparo para este teste. A Figura 16, do lado esquerdo, ilustra a localização desta linha de disparo, que contém 342 disparos. Construiu-se um TTI RTM com a equação de onda quase-P pura proposta e migraram-se estes disparos. A Figura 16, do lado direito, mostra o resultado de imagem usando o algoritmo elíptica aqui descrito, que se sobrepõe com o modelo de densidade. Está claro que a imagem gerada com a nova equação de onda quase-P corresponde ao modelo de densidade muito bem e apresenta resultado limpo e preciso.

[0068] Método numérico

[0069] A Equação 7 é resolvida numericamente para obter componentes de campo de onda para cada localização predefinida dentro da subsuperfície em diferentes tempos de propagação de onda. Isto é feito por um processo de estimativa numérica. O processo numérico para resolver a Equação 7 é relativamente simples de implementar. A solução numérica das Equações 10 e 11 é similar, mas o valor dos escalares será diferente porque o operador diferencial é diferente quando comparado com a Equação 7. O processo tem as etapas S1 a S3 para determinar o campo de onda, tal como estabelecido abaixo e ilustrado na Figura 1: 51. Computar o gradiente do campo de onda atual. O campo de onda atual é, por exemplo, o campo de onda determinado para um intervalo de tempo anterior. Isto é, inicialmente, o campo de onda de cada tempo limite zero. 52. Computar o operador escalar. Isto é feito como indicado na Equação. 5. O operador S é computado em cada localização espacial utilizando o campo de gradiente da etapa S1 para obter as componentes de direção de fase (ignorando o sinal após a projeção de coordenadas locais). Multiplicar o operador escalar e o quadrado da velocidade com o gradiente de campo de onda. Estes são todos valores numéricos. 53. Computar a divergência do resultado da etapa 2, para determinar o campo de onda no dado instante de tempo. Um valor da taxa de variação no tempo a partir do campo de onda atual é obtido, o qual, por sua vez, é usado para obter o campo no novo tempo tipicamente usando um método integral.

[0070] O campo de onda determinado da etapa S3 é usado como o campo de onda atual na etapa S1 do seguinte intervalo de tempo de propagação, como indicado pelo ciclo na Figura 1. As etapas S1 a S3 são repetidas em etapas sucessivas de tempo de propagação, de modo que um campo de onda preciso com relação a espaço e tempo pode ser determinado.

[0071] Como visto na Figura 1, uma etapa inicial S0 do processo pode estar presente para o fornecimento do campo de onda inicial. Uma outra etapa S4 é também mostrada na Figura 1, em que o campo de onda determinado é usado para migrar dados sísmicos. Isto é feito uma vez que o cálculo do campo de onda em cada localização espacial para todos os intervalos de tempo foi concluído. Algoritmos de migração padrão são dispostos de modo a usar o campo de onda ou as suas componentes. O campo de onda apropriado para um tempo de viagem desejado e a localização de interesse na subsuperfície podem ser obtidos a partir do método de solução acima.

[0072] Da mesma forma, a Equação 8 pode ser resolvida por um processo de estimativa numérica com as etapas T1 e T2 como segue: T1. Computar o gradiente do campo de onda atual. O campo de onda atual é, por exemplo, o campo de onda determinado para um intervalo de tempo anterior. Isto é, inicialmente, o campo de onda no limite no tempo zero. T2. Computar o operador escalar S como visto na Equação 5 para cada localização espacial utilizando o campo de gradiente da etapa T1 para obter as componentes de direção de fase (ignorando o sinal após a projeção de coordenadas locais). Multiplicar o operador escalar e o quadrado da velocidade com o gradiente de Laplace do campo de onda. Um valor da taxa de variação no tempo é obtido, o qual, por sua vez, é usado para obter o campo no novo tempo tipicamente usando um método integral.

[0073] Os cálculos de gradiente dentro das etapas relevantes S1 a S3 ou T1 a T2 dos métodos acima podem ser realizados utilizando algoritmos de diferenças finitas padrão, ou alternativamente, usando métodos de Transformada Rápida de Fourier (FFT) padrão.

[0074] A determinação do campo de onda, quando usando a técnica numérica de diferenças finitas custa, no máximo, o dobro do que usando a solução de equação de onda acústica isotrópica padrão. Isto é muito mais rápido do que outras abordagens de propagação da onda anisotrópica.

[0075] A fim de obter o gradiente e a divergência nas etapas S1 a S3 acima, as primeiras derivadas (del e divergência) do campo de onda precisam ser calculadas, e um esquema numérico otimizado pode ser utilizado para calculá-las simultaneamente, o que pode ser eficaz. Da mesma forma, a fim de obter o Laplaciano (del quadrado) na Equação 8, são necessárias as primeiras e segundas derivadas e um esquema numérico otimizado poderia ser utilizado onde a primeira e a segunda derivada do campo de onda poderiam ser calculadas simultânea e eficientemente. Por exemplo, pode-se utilizar a transformada rápida de Fourier (FFT) para computar essas derivadas espaciais. Neste caso, precisa-se apenas de uma FFT direta e duas FFT inversas para se obter a primeira derivada e a segunda derivada simultaneamente. Novamente, os algoritmos de FFT padrão são adequados.

[0076] Quando usando a técnica numérica FFT, a determinação do campo de onda é realizada a um custo numérico adicional de apenas 50% a mais do que a solução da equação de onda acústica padrão para um meio isotrópico. Nota-se que o custo numérico aumenta minimamente quando movendo de meio VTI para TTI ou para TOA.

[0077] A Equação 8 utiliza o operador de Laplace diretamente. Portanto, um esquema numérico eficiente poderia ser mais fácil de implementar. Comparando-se a Equação 7, a Equação 8 tem o mesmo comportamento cinemático, mas diferentes efeitos de amplitude.

[0078] O método integral usado para obter o campo de onda nas etapas S3 ou T2 pode ser um método numérico integral de tempo padrão, por exemplo, uma precisão de segunda ordem do esquema de diferenças finitas de integração temporal ou método de expansão rápida (REM) (Kosloff e outros, 1989).

[0079] Em resumo, a equação de onda (Equação 8) tem a forma —= v v2SV2u.

Numericamente, o campo de onda depende do ponto espacial X e do tempo t, como poderia ser expresso como U(x, t). A tarefa de resolver o campo de onda envolve o uso do campo de onda em uma amostra de tempo atual para computar o campo de onda da próxima amostra de tempo. Assim, para o lado esquerdo da equação (segunda derivada no tempo), o método usado pode ser um esquema padrão referido como o método Integral. Para o cálculo desta integral, precisa-se saber o valor do lado direito da equação, que contém derivada espacial dos campos de onda (primeira ordem na Equação 7, e segunda na Equação 8). Os esquemas numéricos convencionais podem ser utilizados para este fim, que é Diferenças Finitas (FD) ou FFT (FD é mais eficiente, e FFT é mais precisa), enquanto nota-se que, em geral, tanto FD quanto FFT são padrão e são algoritmos bem compreendidos.

[0080] Voltando agora à Figura 2, é representado um dispositivo de computador para executar os métodos de determinação dos campos de onda ou suas componentes, conforme apresentado acima. O dispositivo de computador 10 tem um dispositivo In/Out dispositivo 11, um microprocessador 12 e uma memória 13. Os programas de computador 14a e 14b são armazenados na memória 13. O programa de computador de campo de onda 14a tem instruções para executar os métodos numéricos para resolver as Equações 7 e 8 para obter o campo de onda em diferentes tempos de propagação e localizações espaciais. O microprocessador 12 está disposto para ler e executar as instruções contidas no programa de computador de campo de onda para determinar o campo de onda. O campo de onda calculada é preferencialmente também armazenado na memória, e é transmitido como entrada para um programa de migração também armazenado na memória, também executável pelo processador 12 de modo a efetuar a migração de dados sísmicos usando o campo de onda calculado. O dispositivo também pode ter uma tela para visualizar os dados armazenados na memória e/ou calculados através dos programas. A migração pode ser realizada como descrito acima neste documento, incluindo a seção de fundamentos.

[0081] O dispositivo In/Out 11 é utilizado para ler e emitir dados a partir do dispositivo de computador. Em particular, os dados sísmicos podem ser recebidos através do dispositivo In/Out, como pode ser obtido em um levantamento sísmico, e os dados armazenados na memória 13.

[0082] Na Figura 3, o aparelho é representado compreendendo um aparelho de levantamento sísmico 1. O aparelho inclui um navio de levantamento sísmico rebocando uma fonte sísmica 4 e detetores sísmicos 5 através de um corpo de água. A fonte sísmica é utilizada para transmitir uma onda sísmica através da subsuperfície 2. A onda interage com uma interface 3 e uma parte da energia é refletida de volta para os detectores. Os detectores são dispostos de modo a detectar a energia recebida no detector. Os detectores podem ser tipicamente usados para obter dados compreendendo registros de amplitudes com relação ao tempo de viagem relativo ao evento de origem que gera a onda sísmica. Eventos de grande amplitude podem, então, ser associados com reflexões a partir de interfaces na subsuperfície. Os dados do levantamento podem ser lidos por um dispositivo de computador, e processados para fornecer uma imagem da subsuperfície para ajudar a revelar a estrutura geológica. Por exemplo, o aparelho pode compreender um dispositivo de computador descrito na Figura 2, e os dados do levantamento podem ser preparados por esse dispositivo computador e processados utilizando o processador. Os dados podem então ser migrados para fornecer relatórios de dados de levantamento sísmico migrados ou imagens, por exemplo, seções sísmicas migradas. Tal sistema pode ser aplicado para fornecer dados como descrito neste documento acima.

[0083] O dispositivo de computador pode ser um dispositivo distribuído em que qualquer um ou mais do dispositivo In/Out 11, microprocessador 12, memória 13 e tela 15 podem ser distribuídos através de diferentes localizações. A comunicação entre eles pode acontecer como indicado, através de uma rede, por exemplo, uma rede sem fio. Os programas, dados de campo de onda e/ou dados migrados podem, em certas modalidades, ser armazenados no meio de armazenamento removível, tal como um cartão de memória ou disco compacto, executável pelo dispositivo de computador e/ou processador mediante a conexão com ele. Um sinal pode ser fornecido, que é comunicado através da rede, contendo os programas, instruções legíveis por máquina dos mesmos, os dados do campo de onda e/ou dados migrados produzidos como descrito acima.

[0084] As soluções de campo de onda a partir das Equações de onda presentemente propostas 7 e 8 podem ser usadas em outras aplicações que exigem o uso de um campo de onda previsto de um único modo de onda, por exemplo, modelagem sísmica, inversão completa da forma de onda, etc.

[0085] Exemplos e resultados

[0086] Vários testes dos modelos de equação de onda das Equações 7 e 8 descritas acima foram realizados, resultados que são vistos nas Figuras 4 a 12. Em todos os seguintes exemplos, aplica-se uma técnica pseudo-espectral para derivadas espaciais e esquema de diferenças finitas com uma precisão de segunda ordem para a derivada temporal. Isto é, FFT é usada para as derivadas espaciais (devido à maior precisão), enquanto um método integral com um esquema de diferenças finitas de segunda ordem é utilizado para a derivada no tempo.

[0087] O tamanho da grade numérica é escolhido para evitar a dispersão espacial, isto é, até ao número de onda de Nyquist, e etapa de tempo de propagação é escolhida para satisfazer a condição de estabilidade.

[0088] O primeiro exemplo é um modelo VTI simples, que é homogêneo com a velocidade vertical definida como 2000 m/s e parâmetros anisotrópicos de Thomsen (1986), como ε = 0,2 e δ = 0,1. A ondeleta de origem para este exemplo é a ondeleta de Ricker com uma frequência máxima de 24 Hz. A grade de computação é um cubo 3D com comprimento de 6,0 km com a amostragem espacial como 15 m em todas as 3 direções. A fonte pontual é colocada no meio da grade. A Figura 4 ilustra uma fatia sísmica 2D, localizada no meio da direção Y, do instantâneo de campo de onda 3D no tempo t = 0,8 s. Obviamente, apenas a onda de pseudo-P aparece e não existe nenhuma onda de cisalhamento. A Figura 5 mostra a imagem correspondente do operador S no mesmo tempo instantâneo. Este operador desempenha um papel crucial neste algoritmo. Observa-se que a imagem não é tão lisa e limpa como o campo de onda. Isto ocorre porque a direção de propagação representada por Vu1 . . . > ..... . R perde a sua precisão em torno das localizações onde a função de gradiente tende para zero. Felizmente, esta precisão perdida pode ser trazida de volta pelo operador composto S| V U, que é a magnitude do termo dentro do operador de divergência antes da aplicação do modelo de velocidade, como definido na Equação 7. O efeito combinado do produto de duas quantidades S| V u| é mostrado na Figura 6. A imagem lisa e limpa demonstra que a imprecisão da direção introduz erros pouco perceptíveis ao cálculo do campo de onda.

[0089] Para validar essa abordagem para um caso de TTI, utiliza-se os mesmos parâmetros de computação que os do modelo VTI no primeiro exemplo e estende-se os mesmos para um modelo TTI simples. Portanto, além dos mesmos parâmetros anisotrópicos VTI, introduz-se também o ângulo de inclinação como 45° e o ângulo de azimute como 0°, ou seja, o eixo de simetria anisotrópica é inclinado, a inclinação definida pelos ângulos de inclinação e de azimute (No caso de VTI, o eixo de simetria é vertical (não inclinado)). A Figura 7 mostra a fatia 2D de campos de onda instantâneos na mesma localização e a fatia de tempo como o primeiro exemplo. Nota- se as rotações da frente de onda que é o resultado de parâmetros TTI. De modo similar, a Figura 8 é a saída correspondente do operador S . A Figura 9 ilustra os efeitos combinados dos produtos dos dois operadores S| V u\.

[0090] Também testou-se esta abordagem em um modelo de sal SEAM TTI. As dimensões do modelo são conhecidas como sendo: nx = ny = 864 e nz = 768. A taxa de amostragem de grade é de 10 m em todas as três direções. Coloca-se a localização de origem na posição (x, y, z) = (17, 23,0, 0,0) km. Usa-se novamente a ondeleta de Ricker como a ondeleta de origem, mas desta vez com uma frequência máxima de 75 Hz. O intervalo de tempo de propagação é de 0,5 ms. A Figura 10 ilustra a linha central do modelo de velocidade 3D; e a Figura 11 mostra a fatia central em linha do instantâneo de tempo 3D em t = 1,6 s, e a Figura 12 mostra o instantâneo de tempo em t = 2,4 s. Nota-se que o propagador manipula os campos de onda complexos muito bem. Ele gera as ondas de transmissão, ondas de reflexão, e ondas refratadas, mas não ondas de cisalhamento.

[0091] Vantagens

[0092] Comparando com as equações de onda tradicionais, esta abordagem tem algumas vantagens óbvias. Primeiro é a simplicidade das Equações 7 e 8. Elas mantêm a mesma forma para VTI, TTI e TOA. Em segundo lugar é a eficiência numérica. Comparando com a propagação de ondas em meio TTI, os esquemas convencionais geralmente exigem 3-5 vezes mais recursos de computador do que para simulações de ondas isotrópicas enquanto este esquema proposto somente introduz 50% de custo adicional. Este é um esquema muito mais eficiente do que os existentes. Além disso, o seu desempenho numérico é quase o mesmo para isotropia transversal, isotropia transversal inclinada, anisotropia ortorrômbica, ou anisotropia ortorrômbica inclinada. Em terceiro lugar, é a estabilidade da equação. Similar aos casos acústicos, as instabilidades fracas de TTI no sistema de equações diferenciais de segunda ordem 2 x 2 convencional não aparecem nesta nova equação. Esta solução é numericamente estável para modelos muito complicados, por exemplo, modelo com estruturas de sal complexas e estruturas de excesso de empuxo com mudanças bruscas de eixo de simetria anisotrópica. Uma vez que apenas uma equação diferencial é usada, o novo esquema proposto é mais eficiente do que o de algoritmos convencionais.

[0093] Em resumo, a presente abordagem para determinar o campo de onda sísmica em uma subsuperfície anisotrópica é: • Mais eficiente • Mais preciso • Sem ruído de cisalhamento • Simples de implementar • Numericamente estável

[0094] Embora tendo-se somente discutido o algoritmo para a onda quase-P, a abordagem proposta pode ser facilmente generalizada para resolver problema propagação de onda quase-SV ou até mesmo problemas de atenuação elástica, etc. Isto pode ter valor no futuro para geração de imagens de onda-S ou geração de onda convertida.

[0095] Esta presente solução fornece uma maneira de resolver as equações pseudo-diferenciais. A equação de atenuação é um exemplo de outra tal equação.

[0096] A abreviatura “eq.” é usada para significar “equação” e os termos são aqui utilizados de forma intercambiável.

[0097] Várias modificações e aprimoramentos podem ser feitos sem abandonar o escopo da invenção aqui descrita. REFERÊNCIAS

[0098] Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media. Geophysics 63, 623 a 631.

[0099] Alkhalifah T. 2000. An acoustic wave equation for anisotropicmedia. Geophysics 65, 1239 a 1250.

[0100] Bakker, P. M., e E. Duveneck, 2011, Stability analysis for acoustic wave propagation in tilted TI media by finite differences: Geophysical Journal International, 185, no. 2, 911 a 921.

[0101] Bube, K. P., T. Nemeth, J. P. Stefani, R. Ergas, W. Liu, K. T. Nihei, e L. Zhang, 2012, On the instability in second-order systems for acoustic VTI and TTI media: Geophysics, 77, no. 5 P. T171 a T186 .

[0102] Cheng J. e S. Fomel, 2013, Fast algorithms of elastic wave mode separation and vetor decomposition using low-rank approximation for transverse isotropic media: Expanded Abstracts of 83rd Ann. Internat. Mtg., Soc. Expl. Geophys.

[0103] Chu C., B.K. Macy e P.D. Anno, 2013, Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method: Geophysical Prospecting: 61(3), páginas 556 a 567.

[0104] Duveneck E. e Bakker P.M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics 76, S65 a S75.

[0105] Etgen J.T. e Brandsberg-Dahl S. 2009. The pseudo-analytical method: Application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation. 79th Annual International Meeting, SEG, Expanded Abstracts, 2552 a 2556.

[0106] Fomel S., Ying L. e Song X., 2012, Seismic wave extrapolation using lowrank symbol approximation: Geophysical Prospecting, Volume 61, Issue 3, páginas 526 a 536.

[0107] Fowler P.J., Du X. e Fletcher R.P. 2010a. Coupled equations for reverse time migration in transversely isotropic media. Geophysics, 75, S11 a S22.

[0108] Fowler, P.J. e R. King, 2011, Modeling and reverse time migration of orthorhombic pseudo-acoustic P-waves: SEG, Expanded Abstracts, 30, no. 1, 190 a 195.

[0109] Guan, H., E. Dussaud, B. Denel, P. Williamson, 2011, Techniques for an efficient implementation of RTM in TTI media: 81st Annual International Meeting, SEG Expanded Abstracts, 3393 a 3397.

[0110] Jin, S., R. Madariaga, J. Virieux e G. Lambaré, 1992, Two-dimensional nonlinear inversion of seismic waveforms, numerical results: Geophysics, 51, 1387 a 1403.

[0111] Kosloff, D., Filho A.Q., Tessmer E., Behle A., 1989, Numerical solution of the acoustic and elastic wave equations by a new rapid expansion method: Geophysical prospecting;37(4):383 a 394.

[0112] Lambaré, G., J. Virieux, R. Mandariaga, e S. Jin, 1992, Iterative asymptotic inversion in the acoustic approximation: Geophysics, 57, 1138 a 1154.

[0113] Liu F., Morton S. A., Jiang S., Ni L. e Leveille J.P. 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media. 79th Annual InternationalMeeting, SEG, Expanded Abstracts, 2844 a 2848.

[0114] Nemeth, T., Wu, C., e Schuster, G.T., 1999, Least-squares migration of incomplete reflection data: Geophysics, 64, 208 a 221.

[0115] Sun R., G.A. McMechan, H. Hsiao e J. Chow, 2004. Separating P- and S- waves in prestack 3D elastic seismograms using divergence and curl: Geophysics, 69, 286 a 297.

[0116] Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954 a 1966.

[0117] Yan J. e P. Sava, 2009. Elasticwave-mode separation for VTI media. Geophysics 74, WB19 a WB32.

[0118] Zhan G., R. Pestana e P. L. Stoffa, 2013, An efficient hybrid pseudospectral/finite-difference scheme for solving the TTI pure P-wave equation: Geophys. Eng. 10 025004

[0119] Zhang, H. e Y. Zhang, 2008, Reverse time migration in 3D heterogeneous TTI media: 78th Annual International Meeting, SEG, Extended Abstracts, 2196 a 2200.

[0120] Zhang H., G. Zhang e Y. Zhang, 2009. Removing S-wave noise in TTI reverse time migration. 79th Annual International Meeting, SEG, Expanded Abstracts, 2849 a 2853.

[0121] Zhang Y., H. Zhang e G. Zhang, 2011. A stable TTI reverse time migration and its implementation. Geophysics 76, WA3 a WA11.

[0122] Zhang, H. e Y. Zhang, 2011, Reverse time migration in vertical and tilted orthorhombic media: SEG, Expanded Abstracts, 30, no. 1, 185 a 189.

[0123] Zhou H., G. Zhang e R. Bloor, 2006a. An anisotropic acoustic wave equation for VTI media. 68th Annual Conference and Exhibition, EAGE, Expanded Abstracts, H033.

[0124] Zhou, H., G. Zhang e R. Bloor, 2006b, An anisotropic acoustic wave equation for modeling and migration in 2D TTI media: 76th Annual International Meeting, SEG, Expanded Abstracts, 194 a 198.

Claims

1. Método implementado por computador de migração de dados sísmicos compreendendo: a) determinar um campo de onda em uma subsuperfície anisotrópica da Terra, caracterizado pelo fato de que compreende: prover uma equação de onda de modo de onda única quase acústica desacoplada com base em parâmetros anisotrópicos espacialmente variados, em que a equação de onda compreende um operador diferencial espacial e um operador escalar, S, e em que todos termos da equação de onda que contém um ou mais dos parâmetros anisotrópicos são separadas do operador diferencial espacial, e resolver numericamente com um processador a equação de onda para determinar o campo de onda na subsuperfície anisotrópica; b) fornecer primeiros dados associados com a subsuperfície anisotrópica da Terra; e c) executar migrações sismicas para processar o primeiro dado utilizando o campo de onda determinado para produzir segundos dados associados com a dita subsuperfície.

2. Método, de acordo com a reivindicação 1, caracterizado pelo fato de que a equação de onda de modo de onda única quase acústica é esfericamente decomposta.

3. Método, de acordo com a reivindicação 2, caracterizado pelo fato de que a equação decomposta compreende um operador diferencial esférico e um operador escalar esférico.

4. Método, de acordo com a reivindicação 1, caracterizado pelo fato de que a equação de onda de modo de onda única quase acústica é elipticamente decomposta.

5. Método, de acordo com a reivindicação 4, caracterizado pelo fato de que a equação decomposta compreende um operador diferencial elíptico e um operador escalar elíptico.

6. Método, de acordo com a reivindicação 1, caracterizado pelo fato de que a equação de onda é

onde

onde u é o campo de onda, t é o tempo, S é o operador escalar, v0 é a velocidade na subsuperfície ao longo do eixo de simetria de anisotropia, nh é a direção de fase horizontal, nz é a direção de fase vertical, e onde ε e δ são parâmetros anisotrópicos.

7. Método, de acordo com a reivindicação 1, caracterizado pelo fato de que a equação de onda é

onde u é o campo de onda, t é o tempo, S é o operador escalar, v0 é a velocidade na subsuperfície ao longo do eixo de simetria de anisotropia, nx e ny são a direção de fase horizontal, nz é a direção de fase vertical, e onde ε e δ são parâmetros anisotrópicos.

8. Método, de acordo com a reivindicação 1, caracterizado pelo fato de que o operador escalar, S, tem componentes baseados em anisotropia de subsuperfície e direção de fase do campo de onda, a direção de fase sendo obtida a partir de um gradiente estimado de um campo de onda predeterminado.

9. Método, de acordo com a reivindicação 8, caracterizado adicionalmente pelo fato de que compreende qualquer um de: - estimar um campo de onda predeterminado; - determinar o gradiente do campo de onda; - fornecer os parâmetros de anisotropia; - usar o gradiente e parâmetros de anisotropia para computar o operador escalar S; - combinar o operador escalar S com o quadrado da velocidade e o gradiente do campo de onda para obter um resultado combinado; e - usando o resultado combinado, determinar um componente do campo de onda.

10. Método, de acordo com a reivindicação 1, caracterizado pelo fato de que a equação de onda tem a forma de

em que u é o campo de onda e t é o tempo.

11. Método, de acordo com a reivindicação 1 caracterizado pelo fato de que a equação de onda é u -V.(v 02 SV u )= 0 01

em que u é o campo de onda, t é o tempo, S é um operador escalar, e v0 é a velocidade na subsuperfície ao longo do eixo de simetria da anisotropia.