ES2869401T3

ES2869401T3 - Análisis de velocidad de migración de ecuación de onda usando deformación de imágenes

Info

Publication number: ES2869401T3
Application number: ES12811440T
Authority: ES
Inventors: Paul Sava; Francesco Perrone; Clara Andreoletti; Nicola Bienati
Original assignee: Colorado School of Mines
Current assignee: Colorado School of Mines
Priority date: 2011-07-12
Filing date: 2012-07-12
Publication date: 2021-10-25
Anticipated expiration: 2032-07-12
Also published as: EP2732312A4; EP2732312B1; US9702997B2; EP2732312A1; PT2732312T; WO2013009944A1; US20140321713A1

Abstract

Un método implementado por ordenador (S2100) de análisis de datos de imágenes sísmicas, comprendiendo el método: recibir un primer conjunto de datos de imagen (S2110), representando el primer conjunto de datos de imagen datos de imagen de una estructura geológica obtenidos de receptores sísmicos en una primera posición relativa a la estructura geológica; recibir un segundo conjunto de datos de imagen (S2110), representando el segundo conjunto de datos de imagen datos de imagen de la estructura geológica obtenidos de receptores sísmicos en una segunda posición relativa a la estructura geológica; y caracterizado por que el método comprende adicionalmente medir un desplazamiento relativo (S2130) de los segundos datos de imagen en comparación con los primeros datos de imagen con respecto a la estructura geológica comparando directamente el primer conjunto de datos de imagen con el segundo conjunto de datos de imagen definiendo una función residual, definiendo una función objetivo usando un operador de penalización y calculando un gradiente de la función objetivo.

Description

DESCRIPCIÓN

Análisis de velocidad de migración de ecuación de onda usando deformación de imágenes

Datos de la solicitud relacionada

Esta solicitud reivindica el beneficio de y la prioridad bajo el 35 U.S.C. §119(e) de la Solicitud de Patente de Estados Unidos N. ° 61/506.968, presentada el 12 de julio de 2011, titulada "WAVE-EQUATION MIGRATION VELOCITY ANALYSIS USING IMAGE WARPING". Adicionalmente, los informes titulados "WAVEFIELD-BASED VELOCITY ANALYSIS BY IMAGE WARPING" (adjunto a la solicitud provisional anteriormente referenciada como el Anexo "A") y "WAVEFIELD TOMOGRAPHY BASED ON LOCAL IMAGE CORRELATION" (adjunto a la solicitud provisional anteriormente relacionada como el Anexo B".

Campo

Un aspecto ilustrativo está relacionado con la formación de imágenes sísmicas por extrapolación de campo de onda.

En algunas realizaciones, se estima un modelo de la Tierra basándose en datos sísmicos de reflexión. La idea general descrita en el presente documento es que, si se obtiene un modelo correcto, entonces la energía reflejada se reubica correctamente en el subsuelo en los lugares donde tienen lugar las discontinuidades de las propiedades del material (es decir, donde existen los reflectores en el subsuelo).

Una idea común usada para esto es el principio de semejanza. Este principio dice que, si el modelo de la Tierra es correcto, entonces hay diferentes experimentos de reflectores de imagen en el mismo lugar. De otra manera, no.

Típicamente, la comparación entre imágenes se hace en recopilaciones de imágenes comunes.

Estas son colecciones de imágenes construidas para muchos experimentos en las mismas posiciones en el subsuelo. Por ejemplo, pueden alinearse 100 como una función de la profundidad (z) en las mismas posiciones horizontales (x,y). Estos 100 experimentos se espera que coloquen los reflectores a la misma profundidad, por lo tanto, las "recopilaciones" son "planas" (todos los experimentos observan la misma Tierra).

Este proceso ignora el hecho de que cuando la velocidad es incorrecta, los reflectores de la imagen se mueven no solo hacia arriba y abajo, sino también hacia los lados. Este movimiento no puede observarse puesto que el análisis se hace en una posición horizontal fijada (x,y). Este proceso también requiere las construcciones de imágenes de un gran número de experimentos para analizar cómo cambian las imágenes con el experimento (es decir, evaluar el "movimiento de recopilación").

El proceso ilustrativo propuesto en el presente documento es similar al proceso anteriormente descrito en el que se crean imágenes para diferentes experimentos y puede esperarse que se observen las imágenes de los reflectores en las mismas posiciones en el subsuelo. En contraste a la posición convencional, no se construyen recopilaciones, sino pares (grupos pequeños) de imágenes que se comparan en todas las posiciones. Por lo tanto, puede aprovecharse el principio de semejanza, pero sin recopilaciones. De esta manera, el procedimiento se vuelve más económico (no es necesario tener todos los experimentos antes de que se realice una evaluación de la precisión del modelo), y también más robusto (pueden usarse todas las imágenes en todas las posiciones, a diferencia de las recopilaciones construidas en algunas posiciones). Puede crearse de manera gradual un modelo mientras se acumulan imágenes para diferentes experimentos; mientras que, con las recopilaciones, todas las imágenes son necesarias antes de que se evalúe el modelo para la imprecisión.

El documento US 6 574 563 B1 desvela la comparación de dos conjuntos (en un lapso de tiempo) de datos de imágenes sísmicas.

El procedimiento de actualización de modelo ilustrativo descrito en el presente documento emplea la definición de una función objetivo que minimiza cuando el modelo es correcto. Se usa un procedimiento numérico para minimizar esta función y en el proceso de actualización del modelo.

Para la comparación de imágenes construidas para diferentes experimentos (vecinos) hay diferentes posibilidades. Una puede ser simplemente correlacionar las imágenes en cada punto en el espacio. Si la velocidad es correcta, entonces la correlación de imagen es máxima en algún lugar a lo largo del reflector analizado. De otra manera, tiene lugar la correlación cruzada máxima lejos del reflector analizado. Puede definirse la función objetivo directamente en esta función de correlación cruzada.

Como alternativa, puede evaluarse un campo de vector (el campo de deformación) que apunta en la dirección de la correlación cruzada máxima. Puede medirse también en la imagen un campo de inmersión, que es simplemente ortogonal a los reflectores. Si la correlación cruzada máxima tiene lugar a lo largo del reflector, entonces el campo de deformación es paralelo con el reflector.

Por lo tanto, los campos de inmersión y de deformación son ortogonales. Puede definirse también una función objetivo en estos dos campos: si los vectores son ortogonales (es decir, la velocidad es correcta), entonces su producto escalar es cero, por lo tanto, la función objetivo es mínima. Como alternativa, si los campos de vector no son ortogonales (es decir, la velocidad es incorrecta), entonces su producto escalar no es cero, por lo tanto, la función objetivo no es mínima. El modelo puede a continuación actualizarse en la dirección que minimiza la función objetivo.

Esta es la esencia del método. La porción realmente original es la función objetivo definida en imágenes cercanas, ya sea directamente a través de correlación cruzada o a través de los campos de vector definidos en la correlación cruzada. Pueden usarse también otras medidas de similitud de imagen, como desconvolución de las imágenes, en lugar de su correlación cruzada.

Antecedentes

La formación de imágenes sísmicas implica la estimación de las velocidades de propagación de onda en el subsuelo a partir de datos sísmicos registrados en la superficie. Las velocidades sísmicas están relacionadas con otros parámetros físicos (por ejemplo, densidad y compresibilidad, que caracterizan la litología de la Tierra), y parámetros de roca-mecánicos (por ejemplo, porosidad y sobrepresión de fluidos, que son cruciales en la ingeniería de yacimientos) (Carcione, 2007).

La formación de imágenes sísmicas incluye la estimación tanto de la posición de las estructuras que generan los datos registrados en la superficie como de un modelo que describe la propagación de las ondas en el subsuelo. Los dos problemas están estrechamente relacionados puesto que es necesario un modelo para inferir la posición de los reflectores. Las ondas registradas en la superficie se extrapolan en un modelo del subsuelo (resolviendo una ecuación de onda) y se correlacionan de manera cruzada con un campo de onda de fuente sintética simulado en el mismo modelo (Claerbout, 1985). Bajo la aproximación de dispersión sencilla, los reflectores están ubicados donde los campos de onda de la fuente y receptor coinciden en el tiempo y en el espacio.

La tomografía de ecuación de onda (T arantola, 1984; Woodward, 1992; Biondi y Sava, 1999) es una familia de técnicas que estiman los parámetros de modelo de velocidad a partir de las señales de ancho de banda finitas registradas en la superficie. La inversión se formula normalmente como un problema de optimización donde el modelo de velocidad correcto minimiza una función objetivo que mide la inconsistencia entre un modelo de prueba y las observaciones. La función objetivo puede definirse en el espacio de datos (inversión de forma de onda completa) o en el espacio de imagen (análisis de velocidad de migración).

Inversión de forma de onda completa (FWI) (Tarantola, 1984; Pratt, 1999; Sirgue y Pratt, 2004) trata el problema de la estimación en el espacio de datos y mide el desajuste entre las observaciones y los datos simulados. La inversión de forma de onda completa tiene como objetivo reconstruir el modelo exacto que genera los datos registrados. Adaptando tanto los tiempos de recorrido como las amplitudes, la inversión de forma de onda completa permite que se consiga alta resolución (Sirgue et al., 2010). No obstante, es necesaria una estimación de la fuente, debe modelarse correctamente la física de la propagación de ondas (por ejemplo, isotrópica frente a anisotrópica, acústica frente a elástica, etc.) y una buena parametrización (por ejemplo, impedancia frente a contrastes de velocidad) es crucial (Kelly et al., 2010). Además, un modelo inicial preciso es clave para evitar saltos de ciclo y converger al mínimo global (en lugar de un mínimo local) de la función objetivo.

Debido a la no linealidad de los campos de onda con respecto al modelo de velocidad, la función objetivo en el dominio de datos es altamente multimodal (Santosa y Symes, 1989), y los métodos de optimización local pueden quedarse atrapados fácilmente en los mínimos locales y fallar al converger al modelo correcto. Esto se cumple particularmente para la inversión de forma de onda completa de reflexión. La inversión de forma de onda completa de refracción se centra en las ondas de inmersión y silencia los datos, manteniendo únicamente la energía de la onda de inmersión (Pratt, 1999). Esto conduce a una función objetivo con mejor comportamiento, pero requiere compensaciones muy largas para registrar la energía refractada. Además, este enfoque limita la profundidad en la que puede esperarse un resultado de inversión robusto.

Migration velocity analysis (MVA) (Fowler, 1985; Faye y Jeannot, 1986; Al-Yahya, 1989; Chavent y Jacewitz, 1995; Biondi y Sava, 1999; Sava et al., 2005; Albertin et al., 2006) define la función objetivo en el espacio de imagen y está basado en el principio de semejanza (Al-Yahya, 1989). Si el modelo de velocidad es correcto, las imágenes desde diferentes experimentos deben ser consistentes entre sí puesto que un único modelo de la Tierra genera los datos registrados. Se calcula normalmente una medida de la consistencia a través de semejanza convencional (Taner y Koehler, 1969) o semejanza diferencial (Symes y Carazzone, 1991). Estas dos funcionales analizan un conjunto de imágenes migradas en ubicaciones fijas en el espacio; consideran todas las tomas que iluminan los puntos bajo investigación. El análisis de velocidad de migración conduce a funciones objetivo suaves y a problemas de optimización con buen comportamiento (Symes, 1991; Symes y Carazzone, 1991) y es menos sensible que la inversión de forma de onda completa al modelo inicial. Por otra parte, puesto que una realización ilustrativa no usa amplitudes en la etapa de formación de imágenes, el modelo estimado tiene una resolución inferior que el resultado de la inversión de forma de onda completa ideal (Uwe Albertin, comunicación personal).

El análisis de velocidad de migración mide la invarianza de las imágenes migradas en una dimensión auxiliar (ángulo de reflexión, toma, etc.) (Al-Yahya, 1989; Rickett y Sava, 2002; Sava y Fomel, 2003; Xie y Yang, 2008) o se centra en un espacio extendido (Rickett y Sava, 2002; Symes, 2008; Sava y Vasconcelos, 2009; Yang y Sava, 2011). Todos estos enfoques requieren la migración de toda la prospección para analizar la curva de movimiento en recolecciones de imágenes comunes o la medición del enfoque en una ubicación espacial específica. La dimensionalidad del espacio de imagen (ampliado) y la complejidad computacional de la etapa de análisis de velocidad aprovechan rápidamente escenarios de caso realistas. Además, debido al requisito de memoria elevada para almacenar la información parcial de cada experimento, únicamente puede considerarse un subconjunto de los puntos de imagen en la evaluación de la función objetivo. Los orificios de iluminación y/o las geometrías de obtención irregulares pueden impactar también la calidad de las recopilaciones de imágenes comunes, pero no se informa estudio sistemático de este problema en la bibliografía para nuestro conocimiento.

Se usa de manera rutinaria la migración de tiempo inverso (Baysal et al., 1983; McMechan, 1983) en la geofísica de exploración debido a su capacidad para manejar correctamente la complejidad total del fenómeno de la propagación de ondas (bajo la suposición de que la física de la propagación de ondas en el subsuelo está modelada correctamente). No obstante, su coste computacional es aún prohibitivo y limita su integración en un bucle de análisis de velocidad de migración que requiere la extracción de recopilaciones de imágenes comunes para análisis de movimiento. Un procedimiento de análisis de velocidad de migración basado en un motor de propagación de onda completa permite el aprovechamiento del fenómeno de onda más compleja (por ejemplo, reflexión de inversión, ondas prismáticas y múltiples) y aumenta la cantidad de información en los datos que puede usarse (Farmer et al., 2006). No obstante, debido al coste intrínseco de la extrapolación de ondas en el dominio del tiempo, debe considerarse un enfoque diferente para una implementación efectiva y eficaz de la creación del modelo de velocidad.

Hasta ahora, se ha usado únicamente la herramienta potente ofrecida por la migración de tiempo inverso para reconstruir una imagen de los contrastes de impedancia en el subsuelo. El modelo de velocidad usado para generar una imagen migrada se obtiene con otras técnicas, por ejemplo, a través de tomografía de tiempo de recorrido (Bishop et al., 1985) o a través de análisis de velocidad de migración de ecuación de onda (Biondi y Sava, 1999; Sava y Biondi, 2004; Shen y Symes, 2008), que está basada en aproximaciones asintóticas y/o en la linealización del operador de ecuación de onda. Estas técnicas no aprovechan toda la información codificada en los campos de onda registrados puesto que son imprecisos en los modelos de velocidad complejos (tomografía de tiempo de recorrido) o no pueden modelar de manera apropiada parte del fenómeno de propagación (métodos unidireccionales). Diseñar un procedimiento de creación de modelo de velocidad basándose en un motor de propagación de onda bidireccional permitirá que una realización ilustrativa aproveche la complejidad completa de los campos de onda y hará la etapa de análisis de velocidad conceptualmente consistente con el algoritmo de formación de imágenes usado.

Debido a que la migración de tiempo inverso existe tanto desde un punto de vista computacional como de almacenamiento, una realización ilustrativa propone analizar la semejanza de grupos pequeños de imágenes migradas de experimentos cercanos, es decir tomas vecinas u ondas de plano con parámetro-rayo similar o ángulo de despegue, en cada punto de imagen iluminado por los experimentos. La migración de tiempo inverso funciona únicamente en configuraciones de toma común y receptor común, por lo que puede preferirse la elección del dominio de toma para realizar el análisis de velocidad de migración con un motor bidireccional. En este punto, usamos la palabra "toma" en su sentido más amplio para incluir recopilaciones de toma sintéticas como fuentes de onda plana (Whitmore, 1995; Liu et al., 2006; Stoffa et al., 2006; Zhang et al., 2005), fuentes codificadas por toma aleatoria (Morton y Ober, 1998; Romero et al., 2000) o cualquier otra fuente de fase/amplitud codificada (Soubaras, 2006; Perrone y Sava, 2011) en el intervalo de migración de perfil de toma. Otras técnicas para la creación del modelo de velocidad en el domino de toma usan diferentes estrategias. Sin embargo, todas ellas consideran la prospección entera y migran todos los datos iniciando el bucle de análisis de velocidad de migración (deVries y Berkhout, 1984; Yilmaz y Chambers, 1984; Al-Yahya, 1989; Chavent y Jacewitz, 1995; Sava y Biondi, 2004; Shen y Symes, 2008; Symes, 2008; Xie y Yang, 2008). Una realización ilustrativa desvela una función objetivo que evalúa el grado de semejanza entre imágenes a través de correlaciones locales en el espacio de imagen y no necesita recopilaciones de imágenes comunes (CIG). Se usa una relación morfológica entre las imágenes de experimentos cercanos para definir una función objetivo que mide desplazamientos en el espacio de imagen.

La metodología descrita en este trabajo sigue el operador de MVA de ecuación de onda linealizada propuesto por (Perrone et al., 2012) pero elimina la linealización del operador de onda con respecto a los parámetros de modelo y, por lo tanto, tiene en cuenta la naturaleza bidimensional de los campos de onda. Una realización ilustrativa calcula el gradiente de la función objetivo con el método de estado adjunto (Plessix, 2006). Una ventaja ilustrativa de la técnica desvelada en el presente documento es que reduce los requisitos de memoria y evita la necesidad de captar movimiento en las recopilaciones, y permite la inclusión de todos los puntos iluminados por los experimentos sísmicos en la creación del modelo de velocidad.

Sumario

La estimación de un modelo de velocidad a partir de datos sísmicos es una etapa crucial para obtener una imagen de alta calidad del subsuelo. La estimación de velocidad normalmente se formula como un problema de optimización donde una función objetivo mide el desajuste entre campos de onda sintéticos y registrados y se usa su gradiente para actualizar el modelo. La función objetivo puede definirse en el espacio de datos (como en la inversión de forma de onda completa) o en el espacio de imagen (como en el análisis de velocidad de migración). En general, el último conduce a funciones objetivo suaves, que son monomodales en una cuenca más amplia acerca del mínimo global en comparación con las funciones objetivo definidas en el espacio de datos. No obstante, el análisis de velocidad de migración requiere la construcción de recopilaciones de imágenes comunes en ubicaciones espaciales fijas y el submuestreo de la imagen para evaluar la consistencia entre el modelo de velocidad de prueba y los datos observados.

De acuerdo con una realización ilustrativa, se describe una función objetivo que extrae la información de error de velocidad directamente en el dominio de la imagen sin calcular recopilaciones de imágenes comunes. Debido a la dimensionalidad del problema, se usan los métodos basados en gradiente (tales como el algoritmo de conjugadogradiente) en el procedimiento de optimización. Para incluir la complejidad completa del campo de onda en el algoritmo de estimación de velocidad, la técnica considera un operador de onda bidireccional (a diferencia del unidireccional), y el operador de formación de imágenes no está linealizado con respecto a los parámetros de modelo (como en el análisis de velocidad de migración de ecuación de onda linealizada), y calcula el gradiente de la función objetivo usando el método de estado adjunto. Se ilustra a continuación la metodología de estimación de velocidad con unos pocos ejemplos sintéticos.

Para establecer un problema de optimización, se define un residual en cualquiera del espacio de datos, r(t), o el espacio de imagen, r(x), que representa una medida indirecta del error en los parámetros de modelo. Puede a continuación definirse una función objetivo como la energía (normal b) del residual en el dominio apropiado. Se han usado otras normas (por ejemplo, la norma li) para aprovechar el concepto de escasez y robustez de ganancia frente al ruido inevitable en los datos (Claerbout y Muir, 1973; Tarantola, 2005).

En el caso sencillo de problemas lineales, la norma b conduce a problemas de optimización cuadrática para los que puede calcularse fácilmente el gradiente de la función objetivo; en contraste, la norma li conduce a funciones objetivo no cuadráticas que son mucho más difíciles de diferenciar. Esta diferencia fundamental hace a la norma b una elección más popular para establecer un problema inverso en el marco de la optimización. A continuación, se presenta una breve revisión de cómo se define el residual; a continuación, se demuestra una medida ilustrativa del error de velocidad en el espacio de imagen y la función objetivo resultante. El cálculo de estado adjunto del gradiente concluye la sección.

Residuales del espacio de datos

El residual r(t) puede definirse de diferentes maneras: en inversión de forma de onda completa, el residual es una función del tiempo, y se define como la diferencia entre los datos sintéticos y los observados (Tarantola, 1984):

r ( t ) = d sin( t ) - d obs(t).

Este residual captura información acerca del tiempo de recorrido y el error de amplitud debido a un modelo de velocidad incorrecto. Los datos observados establecen una referencia de tiempo y el desajuste en el tiempo de recorrido puede retroproyectarse directamente en el modelo de velocidad: un retardo de tiempo positivo corresponde a un modelo de velocidad que es demasiado rápido; un retardo de tiempo negativo corresponde a un modelo de velocidad que es demasiado lento. Puede obtenerse una función residual diferente calculando la correlación cruzada de los datos sintéticos y observados y penalizando los retardos distintos de cero de la correlación cruzada (van Leeuwen y Mulder, 2008):

r ( t ) = ^ P ( t) ( ( d sm( t) *T d obs( t )).

T

El asterisco indica la operación de correlación cruzada, la variable ^tindica el retardo de correlación cruzada en el dominio del tiempo, el operador de penalización P(^t ) puede ser cualquier función monotónica de ^tque sea cero en ^t= 0. En este caso, el foco se encuentra únicamente en la información de tiempo de recorrido y no en el error de amplitud en los datos. Si el modelo de velocidad es incorrecto, la correlación cruzada de los datos sintéticos y observados se desplaza del retardo cero e indica el signo del error de velocidad. El operador de penalización aniquila la señal alrededor de retardo cero y destaca el desajuste de tiempo de recorrido.

Residuales del espacio de imagen

En el análisis de velocidad de migración, el residual r(x) se calcula usando los campos de onda de fuente y receptor extrapolados en el modelo de subsuelo, y, por lo tanto, se define en el espacio de imagen. Se informan diversas medidas de velocidad en la bibliografía, la mayoría de ellas están basadas en la curvatura del movimiento observado en las recopilaciones de imágenes comunes en diversos dominios (compensación, ángulo de reflexión, etc.). En la actualidad, se ha usado el criterio de enfoque (Shen y Symes, 2008; Symes, 2008; Sava y Vasconcelos, 2009) para definir un residual en el dominio de la imagen que representa un intermediario para el error de velocidad. Tanto la construcción de recopilaciones de imágenes comunes como la medida de enfoque requieren la migración de todos los datos que iluminan la porción del modelo bajo investigación.

Una función objetivo ilustrativa usada con la técnica desvelada no requiere la construcción explícita de recopilaciones de imágenes comunes y opera directamente en las imágenes migradas. El diseño de esta función objetivo continúa de trabajo anterior sobre la ecuación de onda linealizada MVA (Perrone et al., 2012), que usa la relación de deformación entre imágenes para construir una perturbación de imagen que mide el error del modelo y controla la actualización del modelo. Para una realización ilustrativa, se estiman los desplazamientos relativos entre imágenes usando correlaciones locales penalizadas (un operador bilineal) y puede integrar esta medida en un algoritmo de inversión de campo de onda completo usando el método de estado adjunto. Un objetivo es evaluar la corrección del modelo de velocidad de un pequeño subconjunto de imágenes de experimentos cercanos, pero incluyendo todos los puntos en la imagen en la función objetivo y sin aumentar el coste computacional. Un objetivo es intercambiar la iluminación completa (necesaria para medir el enfoque) en una única ubicación espacial para la apertura del subsuelo (es decir, todos los puntos en la imagen iluminados por las tomas bajo análisis). Se revisa e introduce la medida de consistencia basándose en el principio de semejanza en Perrone et al. (2012), que define un residual r(x) en el espacio de imagen, y a continuación formula una función objetivo como la norma I² del residual.

Existen al menos dos opciones para medir la similitud entre dos imágenes. La primera es calcular la diferencia puntual, caso en el que la función residual es

r ( x ) = Ri_l( x ) - R i (x), ( 1)

donde Ri (x) y R⁺¹ (x) son las imágenes migradas para la toma de orden i y (i + 1), respectivamente. Las dos tomas se supone que están cerca y para la imagen de la misma área en el subsuelo. Plessix (2006) usa la energía de la diferencia de imagen como un término de regularización para la inversión de forma de onda completa y no como una función objetivo independiente. Obsérvese que, si el modelo de velocidad usado en la migración es gravemente impreciso, puede hacerse el salto de ciclo en el espacio de imagen o, en otras palabras, la diferencia entre las dos imágenes puede producir dos eventos si la posición del mismo reflector en las dos imágenes cambia más de un cuarto de una longitud de onda. Este problema es análogo a lo que ocurre en la inversión de forma de onda completa en el espacio de datos.

Una solución alternativa es calcular correlaciones locales de las dos imágenes en cada punto x y penalizar de acuerdo con la inmersión local estimada de una de las imágenes. Este enfoque mide el desplazamiento relativo de una imagen con respecto a la otra directamente en el espacio de imagen y puede usarse para análisis de velocidad. Nuestra metodología es la dual en el dominio de la imagen de la estrategia de correlación de datos desarrollada por van Leeuwen y Mulder (2008) en el marco de la inversión de forma de onda completa. El residual se expresa como

donde

es la correlación local de las imágenes Ri (x) y R⁺¹ (x) (Hale, 2007b) y P (x, A) es un operador de penalización que destaca el desplazamiento de la correlación local a lo largo de la dirección de inmersión. El vector A indica el retardo de correlación en el espacio de imagen y el índice i explora la posición de toma. Las ventanas gaussianas locales w (x) centradas alrededor de los puntos de imagen x ponderan las imágenes de entrada y permiten el cálculo de correlaciones cruzadas.

La condición de formación de imágenes

El residual r (x) en la ecuación 2 se calcula a partir de las imágenes migradas Ri (x) y R⁺¹ (x). Se oculta la dependencia de los campos de onda de fuente y receptor en la condición de formación de imágenes que usamos para construir la imagen migrada. Una imagen del subsuelo R(x) se calcula de manera convencional como la correlación de tiempo de retardo cero de los campos de onda de fuente y receptor:

R {x ) = ^us{x,t)ur{x,t)dt, ( 4)

T

donde T es el intervalo de tiempo de registro en los datos, o la correlación de frecuencia de retardo cero de los campos de onda de fuente y receptor:

R(x) = ^us(x,a>)ur (x,ü))dü), ( 5)

a

donde la barra superior indica el operador de conjugación compleja y O indica el ancho de banda de la frecuencia de las señales.

Función objetivo de correlación de imagen

Puede restablecerse el principio de semejanza como sigue: si el modelo de velocidad es correcto, dos imágenes de experimentos cercanos construyen la imagen final a lo largo de la dirección de los reflectores (Perrone et al., 2012). Esto es equivalente a la suposición de la norma que, si el modelo es correcto, el cubo de imagen previamente apilado es invariante con respecto a la posición de la toma, es decir, no se recoge movimiento en la imagen común del dominio de toma (Xie y Yang, 2008).

La inmersión estructural es un atributo que se extrae de manera común en el procesamiento sísmico (van Vliet y Verbeek, 1995; Fomel, 2002; Hale, 2007a). La similitud de dos imágenes a lo largo de la inmersión estructural puede medirse por medio de correlaciones locales (Hale, 2007b). Esta idea puede ilustrarse en un modelo sencillo con una interfaz de densidad horizontal sencilla (véase la Figura 1); el modelo de velocidad correcto es constante y equivale a 200 km/s. La Figura 2 muestra el campo de vector de desplazamiento de inmersión y aparente (flechas blancas y negras, respectivamente) en las ubicaciones espaciales particulares. Obsérvese que, los dos vectores son ortogonales cuando el modelo de velocidad es correcto.

El cálculo del campo de desplazamiento requiere un procedimiento altamente no lineal, que extrae el máximo de las correlaciones locales en cada punto de la imagen (Hale, 2007b). Se reformula el problema de la evaluación de la ortogonalidad entre la inmersión y el desplazamiento usando las correlaciones locales (como se define en la ecuación 3) penalizadas a lo largo de la dirección de inmersión (Figura 3). Si el modelo de velocidad es correcto, el máximo de la correlación radica a lo largo de la pendiente del reflector (Figura 3(a)); de otra manera, se observa una desviación que representa el desplazamiento relativo de una imagen con respecto a la otra (Figura 3(b)). Para medir el desplazamiento entre las dos imágenes, se considera el operador de penalización

donde A es el vector de retardo de correlación y el vector de inmersión v (x) es la normal al reflector en x. El operador de penalización P (x, A) es una función lineal en la dirección de inmersión (normal al reflector) y las líneas de isopenalización son paralelas al reflector. El operador de penalización es idénticamente cero en la dirección tangencial al reflector. La correlación local de las dos imágenes es una función par con respecto a la dirección del reflector: es así debido a que las dos imágenes iluminan la misma área del modelo y por lo tanto son similares. El operador de penalización cambia la simetría de la correlación local: si el modelo de velocidad es correcto, el producto entre P (x, A) y c (x, A) es una función impar en la dirección de inmersión y apilando a través de los retardos de correlación obtenemos cero (Figura 3(a)); por otra parte, si el modelo de velocidad es incorrecto, el máximo de la correlación local se desvía de la dirección del reflector (Figura 3(b)) y se potencia por el operador de penalización. La pila a través de los retardos de correlación no suma cero puesto que la función no es impar en la dirección de inmersión y obtenemos una medida de desplazamiento relativo entre las dos imágenes. Puesto que el máximo de la correlación local es positivo, el signo del operador de penalización en la dirección de inmersión indica la dirección del desplazamiento relativo de una imagen con respecto a la otra. La Figura 4 muestra las correlaciones locales, operadores de penalización y correlaciones locales penalizadas para una localización particular sobre el reflector en estos tres casos en la Figura 2. Las filas de la Figura 4 hacen referencia a diferentes modelos (de arriba a abajo: uno demasiado alto, el correcto y un modelo de velocidad demasiado baja). El punto en la imagen tiene la posición lateral x = 3,5 km y la vertical que depende del reflector representado.

La Figura 4(a), 4(d) y 4(g) muestra las correlaciones locales calculadas de dos tomas cercanas; obsérvese el desplazamiento en el pico de la correlación local y la orientación de los paneles de correlación. Las Figuras 4(b), 4(e) y 4(h) muestran los operadores de penalización que se usaron para destacar la inconsistencia entre la inmersión medida en la imagen y las correlaciones locales. La orientación de los operadores de penalización, de hecho, depende de la inmersión medida. El ahusamiento aplicado reduce los artefactos de truncamiento de los bordes de los paneles de correlación local. Las Figuras 4(c), 4(f) y 4(i) muestran las correlaciones penalizadas. El desplazamiento aparente es ahora claro como un desequilibrio de signo en los paneles. Cuando el modelo de velocidad es correcto, la correlación penalizada es una función impar (es decir, valor medio cero) (Figura 4(f)); por otra parte, los paneles pierden su simetría cuando el modelo es incorrecto (Figuras 4(c) y 4(i)). El signo del valor medio de los paneles penalizados mide el desplazamiento aparente entre las dos imágenes.

La función objetivo se define como:

donde m representa los parámetros de modelo. La Figura 5 muestra los valores de la función objetivo para diferentes perturbaciones constantes del modelo para el ejemplo sencillo mostrado en la Figura 1. La función objetivo es suave y convexa en el intervalo de errores considerado.

Cálculo del gradiente con el método de estado adjunto

Debido a la dimensionalidad del modelo, se trata el problema de la optimización usando un método basado en gradiente local (descenso más pronunciado, gradiente conjugado, etc.) (Vogel, 2002). El cálculo de la derivada de las variables de estado con respecto a los parámetros de modelo (derivadas de Fréchet) no es práctico (o incluso posible) debido al gran número de dimensiones del espacio del modelo. El método de estado adjunto (Lions, 1972; Plessix, 2006) es un algoritmo eficaz que calcula el gradiente de una función objetivo, que depende de algunas variables que describen el estado del sistema físico bajo análisis (variables de estado), sin calcular las derivadas de Fréchet de estas variables con respecto a los parámetros de modelo.

El método de estado adjunto comprende cuatro etapas: el cálculo de las variables de estado, el cálculo de las fuentes adjuntas, el cálculo de las variables de estado adjunto y, finalmente, el cálculo del gradiente de la función objetivo. La primera y tercera etapa (la evaluación de las variables de estado y estado adjunto) requieren la solución de la ecuación de onda que rige la física del problema. El método de estado adjunto es computacionalmente eficaz puesto que únicamente se requieren dos simulaciones de campo de onda en cada iteración para calcular el gradiente de la función objetivo en lugar de una simulación para cada parámetro en el modelo. Las derivadas de Fréchet de las variables de estado con respecto a los parámetros de modelo son objetos extremadamente grandes: se definen en el espacio dado el producto cartesiano entre el espacio de las variables de estado y el espacio del modelo. No son prácticas de calcular y son imposibles de almacenar en memoria para escenarios de caso realistas en geofísica de exploración. No obstante, describen la sensibilidad de cada variable de estado a cambios en cada parámetro de modelo y son una herramienta potente para análisis de resolución.

A continuación, la variable de estado genérica se denomina usando la letra u, las fuentes adjuntas se indican como g, y las variables de estado adjunto se designan por la letra a. Para tomografía de campo de onda, hay dos fuentes adjuntas por experimento i: el campo de onda de fuente us,¡ (x, t) y el campo de onda del receptor ur,i (x,t), donde i representa el índice de la toma.

Las variables de estado u se obtienen de la solución al siguiente problema

L(m)u = / , (6 )

donde f es el vector de término de fuente, L(m) es el operador de modelación directo, y m indica los parámetros de modelo. L puede ser un operador de onda unidireccional o bidireccional; por lo tanto, puede usarse cualquier continuación hacia abajo (Stolt, 1978; Gazdag y Sguazzero, 1984; Stoffa et al., 1990) migración de tiempo inverso (Baysal et al., 1983; McMechan, 1983) en el procedimiento de tomografía de campo de onda. En este punto, consideramos un operador de onda bidireccional (d'Alambert o Helmholtz), y m representa la lentitud al cuadrado. Como se muestra a continuación, la elección de lentitud al cuadrado es conveniente puesto que simplifica el cálculo final del gradiente y nos permite obtener una expresión que es independiente de los parámetros de modelo en la iteración actual.

Las fuentes adjuntas g son las derivadas parciales de la función objetivo con respecto a las variables de estado:

53

(7)

8 ~ du'

La expresión real de las fuentes adjuntas depende de la función objetivo. El diseño particular de la función objetivo impacta las fuentes adjuntas y caracteriza la estrategia de tomografía de campo de onda particular propuesta. En las siguientes secciones, las expresiones para las fuentes adjuntas se derivan para una función objetivo basándose en correlaciones de imagen locales (ecuación 6).

Se resuelve el problema adjunto y se calculan las variables adjuntas usando las fuentes adjuntas g como el término de fuerza:

L*(m)a = g, (8)

donde L* (m) es el operador adjunto a L(m). Finamente, se proporciona el gradiente por el producto interno

03 cL

— = - { a ,— (9 )

cin \ dmi

Por ejemplo, para el operador de Helmholz parametrizado en términos de lentitud al cuadrado, L = -v 2-/?76u2, la derivada dL I dL\ 2

a _{,— ) =}- _ü) (a,u) parcial devuelve un factor de escalamiento sencillo -cu2, y el producto interno \ ' d 'Jm" l l / puede observarse también como la correlación de retardo cero de u y a, similar al procedimiento empleado en FWI convencional.

Fuentes adjuntas

_ 03

Como se describe en la sección anterior, las fuentes adjuntas du son las derivadas de la función objetivo con respecto a las variables de estado. En nuestra función objetivo, se considera el campo de inmersión una función de variables de estado que varía lentamente y desprecia la derivada del operador de penalización con respecto a u. Los detalles de la derivación de las expresiones para las fuentes adjuntas se detallan más adelante; en este punto, el resultado final se informa como:

donde gs,¡ = gs,i (x, q, t) y grj = gr,¡ (x, q, t) indican las fuentes adjuntas para el campo de onda de fuente y receptor de la toma de orden i, respectivamente, y r = r (x) se define como en la ecuación 2.

Las expresiones de las fuentes adjuntas en las ecuaciones 10 y 11 son bastante complejas y necesitan explicarse en algún detalle. Obsérvese también que las fuentes adjuntas dependen de 3 variables (x, n, t). El vector x representa el espacio físico, donde se define el modelo y la imagen; la variable t identifica el eje de tiempo; q es un vector auxiliar definido en el espacio físico que abarca la ventana local alrededor de cada punto de imagen.

La fuente adjunta se obtiene ensanchando el valor del residual en cada punto en una ventana alrededor del punto de la imagen y ponderándola mediante la integral a través de A, que representa una convolución local de la imagen y el operador de penalización. El residual en la ecuación 2 mide el desplazamiento relativo de una imagen con respecto a la otra; la fuente adjunta estima de esta manera la curvatura del movimiento en las recopilaciones de imagen común del dominio de las tomas en cada punto de imagen y escala los campos de onda de fondo mediante este valor. Las ecuaciones 10 y 11 dependen de dos imágenes, R- i (x) y R+i (x), y pueden sufrir de salto de ciclo si la distancia entre las dos tomas y la complejidad del modelo provocan que los reflectores reconstruidos estén más distantes de un cuarto de una longitud de onda. Para superar este problema, se usa la imagen central R¡ (x) para calcular las fuentes adjuntas.

La dependencia del operador de penalización P (x, A) en las variables de estado u pasa a través de la definición del vector de inmersión u (x), que es normal a los reflectores en cada punto en la imagen. El campo de inmersión define el conjunto de frente de onda de la función que representa los eventos de reflexión en la imagen (Chang et al., 1987); puede escribirse el vínculo entre la imagen (y los campos de onda) y el campo de inmersión aprovechando el concepto de tensor cuadrado de gradiente (van Vliet y Verbeek, 1995): el vector de inmersión es el vector propio asociado con los valores propios más grandes asociados. No hay relación lineal sencilla entre el vector de inmersión y los campos de onda, es decir, las variables de estado. Despreciando este término en el cálculo de las fuentes adjuntas puede introducir un error del mismo orden de magnitud del término considerado. Un estudio minucioso de SP/ Su y su impacto en el cálculo tanto de las fuentes adjuntas como del gradiente de la función objetivo es un objeto para investigación futura.

En la siguiente prueba, se consideró la interfaz de densidad sencilla en la Figura 1 y tres tomas 60 m separadas empezando desde x = 4 km. Los datos se calcularon en un medio de velocidad de 2,0 km/s constante usando un esquema de diferencia finita bidireccional. Se migran los datos con tres modelos de lentitud diferentes y se realiza análisis de sensibilidad para la metodología. Se muestran los residuales, la función objetivo como una función de la perturbación de modelo, y el gradiente con respecto a los parámetros de modelo.

Reflector plano en un medio de velocidad constante

La Figura 6 muestra las imágenes migradas de una toma ubicada a x = 4 km usando tres modelos de velocidad diferentes. Obsérvese que, con una única toma, la imagen de un reflector horizontal se curva arriba o abajo dependiendo del signo del error de velocidad. El desplazamiento relativo entre diferentes imágenes se mide por el residual definido en la ecuación 2. La Figura 7 muestra los residuales calculados usando dos imágenes de tomas adyacentes. Obsérvese el cambio en el signo a través del punto de reflexión de compensación cero y cómo cambia el signo con el error de velocidad. El desplazamiento entre las imágenes nos indica que el modelo de velocidad contiene errores, pero no están relacionados de manera inequívoca con el signo del error. Pueden combinarse los residuales y obtener una estimación de la curvatura de las recopilaciones de imagen común del dominio de las tomas (Figura 8) que está más directamente conectada con el signo del error de velocidad. El volteo de signo en los bordes del reflector es debido a la apertura limitada del experimento y al hecho de que las tres imágenes no solapan a través de toda la extensión del reflector.

En la Figura 9, se muestra el gradiente calculado a partir de la función objetivo de correlación. Se consideraron tres tomas para calcular los residuales y para usar únicamente los campos de onda del experimento central para calcular el gradiente. Obsérvese la correlación entre el signo y el patrón de la estimación de la curvatura en la Figura 8 y el gradiente calculado de la función objetivo. El gradiente de la función objetivo representa la dirección del aumento máximo de la misma función y su signo depende de la parametrización elegida. Se calcula a continuación el gradiente con respecto a la lentitud al cuadrado; esto significa que un gradiente positivo apunta hacia modelos más lentos mientras que un gradiente negativo apunta hacia modelos más rápidos.

Ejemplos

Se calculan los gradientes para la función objetivo basándose en correlaciones de imagen locales usando el método de estado adjunto y se muestran unas pocas pruebas de inversión. Se analizan 3 casos de complejidad creciente. En primer lugar, se considera una interfaz de densidad horizontal (Figura 1) en un modelo de lentitud homogénea. Se usan 30 fuentes separadas 200 m de manera equitativa en la superficie y se ubican los receptores en cada punto de cuadrícula a z = 0. Esto demuestra cómo se actualiza el modelo a través de la inversión por medio de un método basado en gradiente y el método de estado adjunto. En segundo lugar, se ejecuta una prueba de inversión en un modelo heterogéneo sencillo con diferentes capas usando 40 tomas separadas 200 m de manera equitativa y los receptores en cada punto de cuadrícula en la superficie. Se usan las recopilaciones de imágenes comunes del dominio de las tomas para evaluar la calidad del resultado y para mostrar que la metodología ilustrativa puede obtener un conjunto de imágenes que es invariante con respecto a la posición del experimento. La prueba final usa el modelo de Marmousi. Los datos se generan con un único código de dispersión y se migran con un esquema de continuación hacia abajo.

Modelo de lentitud homogéneo

Se considera el modelo sencillo en la Figura 1 con una única interfaz de densidad horizontal. Los datos se generan usando un código de diferencia finita del dominio del tiempo con condiciones de límite absorbentes implementadas en todos los lados del modelo. Se añade ruido gaussiano aleatorio a los datos y la relación de señal a ruido es igual a 10. El modelo correcto (Figura 10) es una capa de lentitud homogénea y el modelo inicial es homogéneo, pero con un valor incorrecto de lentitud (Figura 11 (a)). Se modelan 30 tomas en la superficie empezando desde x = 1,01 km; el espaciado entre las tomas es de 200 m y los 400 receptores están ubicados en cada punto de cuadrícula en la superficie.

En la imagen obtenida con el modelo inicial (Figura 11(b)), el reflector está situado incorrectamente debido al valor erróneo de los parámetros de modelo; las recopilaciones de imagen común del dominio de las tomas calculadas para los propósitos de ilustración muestran que diferentes tomas representan el reflector a diferentes profundidades (Figura 11(c)). Se ejecutaron cinco iteraciones de descenso más pronunciado (Vogel, 2002) y se reconstruyó el modelo de lentitud por encima del reflector (Figura 12(a)). La imagen migrada obtenida con el modelo final (Figura 12(b)) se desplaza hacia arriba como resultado del valor de lentitud superior. Las recopilaciones de imágenes comunes en la Figura 12(c) muestran el enfoque obtenido con el modelo reconstruido. Los bordes izquierdo lejano y derecho lejano del reflector muestran la curvatura residual, que es debida a la ausencia de iluminación y de restricciones de los datos. La ausencia de iluminación conduce a un valor superior de lentitud en la parte central del modelo: la inversión compensa la lentitud inferior en los lados con una lentitud superior en el medio. Obsérvese la rápida reducción de la función objetivo a medida que se aplanan los recopiladores (Figura 13).

Modelo heterogéneo lateralmente sintético

Se usa el modelo heterogéneo sintético en la Figura 14 para probar el algoritmo de inversión. El modelo tiene pocas capas con diferentes inmersiones y una estructura inclinada con una inversión en la tendencia de lentitud que de otra manera sería decreciente. Se generan datos acústicos completos (es decir, no de Born) con condiciones de límite absorbente (sin múltiplos de superficie libre) y el ruido gaussiano aleatorio con una relación de señal a ruido igual a 10. Una realización modela 40 tomas separadas 200 m de manera equitativa y los receptores en cada punto de cuadrícula en la superficie. El modelo inicial está altamente suavizado en la Figura 14 y el modelo de lentitud inicial para migración obtenido (Figura 15(a)). La Figura 15(b) muestra la pila de las imágenes migradas obtenidas del modelo inicial; si se compara el modelo correcto en la Figura 14 y la imagen migrada, puede observarse la situación incorrecta de las interfaces de reflexión.

Se construyen recopilaciones de 10 imágenes comunes del dominio de las tomas para evaluar el enfoque global de la imagen. Estas recopilaciones son simplemente la yuxtaposición de las imágenes migradas en posiciones laterales fijadas (de x = 0 km a x = 9 km cada 1 km).

La Figura 15(c) muestra tales recopilaciones para el modelo inicial: la posición de los reflectores cambia como una función del experimento (las imágenes no satisfacen el principio de semejanza) y el modelo, por lo tanto, no es preciso. Después de 15 etapas de tomografía de forma de onda, se recupera el modelo en la Figura 16(a). El modelo es aún bastante suave pero la imagen migrada muestra una mejora perceptible en la posición de los reflectores (Figura 16(b)), y las recopilaciones de imagen común del dominio de las tomas en la Figura 16(c) indican de manera clara mejor enfoque. Además, el signo de los reflectores es ahora consistente con los contrastes de lentitud reales en el modelo exacto. Se impone conformidad entre la actualización del modelo y las capas en la imagen por medio de un operador de suavizado orientado a la estructura, que dirige el gradiente de la función objetivo suavizando a lo largo de la pendiente del reflector. Este enfoque es conceptualmente similar a la inversión escasa propuesta por Ma et al. (2010).

La Figura 17 muestra la evolución de la función objetivo con iteraciones. Obsérvese la reducción monotónica en el valor del residual y el aplanamiento cuando converge el algoritmo. La inversión también muestra una rápida convergencia en las primeras iteraciones de descenso más pronunciado: el suavizado orientado a estructura fuerza el gradiente para que se adapte a la geometría de las capas y acelera la convergencia a un modelo aceptable. A pesar de la simplicidad del modelo, la imprecisión de la cinemática de onda y la ausencia de energía de inversión (ondas de división) hacen a las técnicas del dominio de datos (tales como inversión de forma de onda completa) ineficaces en este caso particular.

Modelo Marmousi

Una prueba de inversión más compleja es configurar usando el modelo Marmousi. En este punto, se modelan los datos usando un código de modelación de dispersión única de diferencia finita (es decir, Born) y condiciones de límite absorbente en los 4 lados del modelo. Se simularon 78 tomas, el espaciado entre las tomas es 0,08 km y la primera fuente se encuentra en x = 0,96 km; se colocaron los receptores en cada punto de cuadricula en z = 0 km. No se manejaron correctamente las inmersiones que entraban en conflicto por una realización ilustrativa puesto que hacen ambigua la definición de los operadores de penalización. Por esta razón, la inversión está restringida a los primeros 1,5 km de profundidad, donde los reflectores son más coherentes.

El modelo inicial es una versión altamente suavizada del modelo de lentitud Marmousi correcto. La Figura 18 muestra el modelo de lentitud inicial, la imagen migrada asociada y las recopilaciones de imágenes comunes del dominio de las tomas extraídas cada 1 km de x = 0 km. Debido a la regularidad de las estructuras geológicas en la parte poco profunda del modelo, la imagen migrada no se ve gravemente distorsionada, no obstante, los reflectores están situados de manera incorrecta y las fallas no están correctamente enfocadas. Las recopilaciones de imagen común del dominio de las tomas (Figura 18(c)) muestran la variación de la profundidad de los reflectores visualizada como una función del experimento, que indica la imprecisión del modelo.

Se ejecutaron 14 iteraciones de descenso más pronunciado y se obtuvieron los resultados en la Figura 19. Se corrige el modelo de lentitud por el procedimiento de inversión y las interfaces visualizadas se mueven hacia la posición correcta. Los planos de falla se enfocan mejor y reducimos los cruzamientos de interfaz en las desconformidades (x = 3 km, x = 4 km y x = 5 km). La recopilación de la imagen común del dominio de las tomas en la Figura 19(b) muestra eventos horizontales más planos que indican un modelo cinemático más preciso y se enfocan eventos adicionales en la parte inferior del modelo (a aproximadamente z = 3 km). Obsérvese también que las recopilaciones de imágenes comunes en partes pobremente o nada iluminadas del modelo contribuyen principalmente al ruido de migración y, por lo tanto, no pueden usarse para evaluar la corrección del modelo de velocidad.

En la Figura 20, la función objetivo muestra una reducción constante. Se calculan los residuales a través de toda la imagen, pero una realización ilustrativa usa únicamente los puntos en los primeros 1,5 km de profundidad para actualizar el modelo de velocidad. Para áreas complejas con inmersiones en conflicto y características geológicas complicadas (donde no hay definición clara del campo de inmersión), un diseño más sofisticado del operador de penalización es clave para obtener residuales significativos y un gradiente fiable. Pueden usarse los valores propios y los vectores propios de los tensores de gradiente cuadrados (van Vliet y Verbeek, 1995) para definir elipses que, su vez, pueden ofrecer un criterio orientado a estructuras para la definición de los operadores de penalización.

Una imagen completa del subsuelo es la superposición de imágenes parciales de experimentos individuales. El principio de semejanza (Al-Yahya, 1989) es un criterio común para evaluar la corrección del modelo de velocidad usado para realizar la formación de imágenes de la prospección: cuando el modelo de velocidad es correcto, todas las tomas ubican los reflectores en la misma posición, es decir, la imagen es invariante a lo largo del eje del experimento. Son necesarias varias tomas para evaluar un error de velocidad en un único punto en el espacio.

Puede evaluarse la invarianza a lo largo del eje del experimento calculando la energía de la primera derivada en esa dimensión. La primera derivada actúa como un operador de penalización destacando y potenciando desviaciones de la dirección horizontal a lo largo del eje de la toma. La pila de la energía de la primera derivada es el operador de semejanza diferencial aplicado directamente en el domino de las tomas (Symes, 1991; Plessix, 2006).

En este punto, se explora una declaración alternativa del principio de semejanza: cuando el modelo de velocidad es correcto, interfieren de manera constructiva las imágenes de diferentes tomas y construyen la imagen perpendicular a la inmersión estructural o, de manera equivalente, paralela a la pendiente del reflector, en cada punto en el espacio de imagen. La inmersión estructural es un atributo extraído de manera común y puede vincularse a la misma imagen por medio de los tensores de gradiente cuadrados (van Vliet y Verbeek, 1995). Desafortunadamente, la relación entre los campos de onda y el campo de inmersión es no lineal: el vector de inmersión representa el vector propio asociado con el valor propio más grande del tensor de gradiente cuadrado. La integración de información acerca de la variación de inmersión con respecto a la perturbación del modelo no es directa, y es necesaria investigación adicional para desarrollar un método eficaz para aprovechar esta información para propósitos de análisis de velocidad. No obstante, puede medirse la semejanza de dos imágenes a través de correlaciones locales penalizadas de manera apropiada de pares de imágenes. Si el modelo de velocidad es correcto, el máximo de la correlación local es a lo largo del reflector en cada punto en la imagen; si el modelo es incorrecto, el máximo se desvía de la pendiente del reflector. El operador de penalización es dependiente del espacio y aniquila los paneles de correlación ortogonales a la inmersión del reflector. Debido a la dependencia del modelo de velocidad, se mide el campo de inmersión en cada iteración tomográfica; la estimación puede llevarse a cabo eficazmente usando tensores de gradiente cuadrados (van Vliet y Verbeek, 1995; Hale, 2007a).

La función objetivo de correlación es superior a la diferencia de imagen en muchos aspectos. En primer lugar, la diferencia entre dos imágenes depende de las amplitudes de la imagen que cambian como una función de la posición de la toma y no puede adaptarse de manera puntual como en la implementación convencional de la inversión de forma de onda completa en el espacio de datos. Los patrones de amplitud afectan al residual en el espacio de imagen y contribuyen de manera eficaz al cálculo de la fuente adjunta, incluso si el modelo de velocidad es correcto y el gradiente de la función objetivo es cero (Mulder y Kroode, 2002). Penalizando correlacionales locales, se reduce la dependencia de la función objetivo en las amplitudes, aumentando de esta manera la robustez y reduciendo la desviación sistemática provocada por las diferencias de amplitud entre las imágenes. Una desventaja del operador de correlación es la pérdida de resolución espacial. La correlación en el dominio espacial es equivalente a la multiplicación en el dominio de la frecuencia/número de onda dual; para señales tanto con soporte espacial finito como con ancho de banda, la multiplicación de los espectros reduce el ancho de banda, es decir, aumenta la anchura de la señal en el dominio espacial. Para aumentar la resolución y precisión de la evaluación de los desplazamientos relativos entre imágenes, la desconvolución es una mejora viable sobre la correlación. Si se supone que dos imágenes de experimentos cercanos están vinculadas por un único desplazamiento espacial sencillo, un enfoque de desconvolución local produciría de manera ideal una función delta espacial limitada en banda que indica la dirección del desplazamiento aparente. Por otra parte, puesto que la desconvolución amplifica el ruido en los datos (debido a la división en el dominio de la frecuencia), es necesario un cuidado adicional para estabilizar el resultado.

Un estudio minucioso del efecto de la distancia de la toma en el gradiente calculado sería beneficioso en relación con el fenómeno, tal como el salto de ciclo, que dificulta muchas estrategias de análisis de velocidad (especialmente en el dominio de datos pero no exclusivamente). El dominio de la imagen es intrínsecamente menos propenso a problemas de salto de ciclo; sin embargo, permanece la cuestión acerca de qué ocurre cuando dos imágenes iluminan áreas del subsuelo bastante diferentes y las correlaciones locales no pueden usarse para estimar un desplazamiento fiable y significativo entre las dos imágenes. A partir de nuestro conocimiento y entendimiento actual, las tomas deben ser lo suficientemente cercanas para proporcionar imágenes del subsuelo comparables y evitar el salto de ciclo en el dominio de imagen. El muestreo de la posición de toma y el patrón de iluminación en el subsuelo pueden crear también un escenario en el que las dos imágenes no solapan en ciertas ventanas locales.

El método ilustrativo está basado en eventos localmente coherentes, tales como reflectores localmente suaves. Las inmersiones, planos de falla y áreas que entran en conflicto donde la definición de un plano de reflexión es ambigua, representan problemas abiertos puesto que en estas áreas no pueden definirse operadores de penalización. A este respecto, la migración de onda de plano puede representar una solución válida debido a la filtración espacial implícita de la imagen: cada onda de plano reconstruye un subconjunto particular de las inmersiones en el modelo, reduciendo por lo tanto la ambigüedad al definir los operadores de penalización.

En contraste con otras técnicas en el espacio de imagen (Bishop et al., 1985; Biondi y Symes, 2004; Lambare et al., 2004; Xie y Yang, 2008; Yang y Sava, 2010, 2011), no se necesitan calcular recopilaciones de imágenes (o imágenes ampliadas). Se proporciona un factor de ahorro de coste adicional mediante la selección de puntos en los reflectores representados (Yang y Sava, 2010) y lejos de áreas complejas (pellizcos o áreas con inmersiones en conflicto). El cálculo de la correlación local se lleva a cabo en cada punto de imagen usando el método eficaz desarrollado por Hale (2006).

La capacidad para extraer información de velocidad de pares de experimentos añade un nuevo grado de libertad para implementar un procedimiento de actualización de modelo. En este punto, la técnica indica dos estrategias posibles. Puede incluir simultáneamente todas las tomas en la prospección en la definición de la función objetivo o puede continuar de manera iterativa y actualizar el modelo de velocidad usando la información obtenida a partir de un único par de experimentos antes de moverse al siguiente par. Las tomas cercanas sondean porciones similares del modelo y proporcionan imágenes comparables; puede usarse la información extraída de un grupo inicial para actualizar el modelo usado para formar imágenes de un segundo grupo de experimentos. Puesto que en cualquier esquema de análisis de velocidad de migración tiene que representarse la prospección entera al menos una vez, la actualización iterativa del modelo a través de las tomas se vuelve rentable si se reduce realmente el número de migraciones globales del conjunto de datos entero.

En 2D, puede definirse fácilmente un orden de los experimentos; en el escenario en 3D general, se tiene un grado extra de libertad. La definición de la función objetivo sigue siendo la misma pero una toma dada puede tener más de 2 experimentos vecinos (la adquisición se define en una cuadrícula bidimensional). Puede generalizarse el concepto de correlación usando la semejanza funcional (Taner y Koehler, 1969), que es no lineal en las señales de entrada y hace el cálculo del gradiente más implicado. Puede analizarse también de manera separada cada par de experimentos dado un modelo de referencia. Es necesaria investigación adicional y pruebas numéricas para evaluar la eficacia y robustez de una estrategia sobre la otra.

Breve descripción de la invención

Los aspectos de la invención se refieren, por lo tanto, a un método implementado por ordenador y/o a un medio de almacenamiento de información legible por ordenador no transitorio para analizar datos de imágenes sísmicas, que comprende:

recibir un primer conjunto de datos de imagen, representando el primer conjunto de datos de imagen datos de imagen de una estructura geológica obtenidos de receptores sísmicos en una primera posición relativa a la estructura geológica;

recibir un segundo conjunto de datos de imagen, representando el segundo conjunto de datos de imagen datos de imagen de la estructura geológica obtenidos de receptores sísmicos en una segunda posición relativa a la estructura geológica; y

medir un desplazamiento relativo de los segundos datos de imagen en comparación con los primeros datos de imagen con respecto a la estructura geológica comparando directamente el primer conjunto de datos de imagen con el segundo conjunto de datos de imagen que define una función residual,

definir una función objetivo usando un operador de penalización, y calcular un gradiente de la función objetivo.

Los aspectos se refieren adicionalmente al elemento estructural de la estructura geológica que comprende una inmersión estructural.

Los aspectos se refieren adicionalmente a la medición de un desplazamiento relativo, usando un ordenador, del primer conjunto de datos de imagen con el segundo conjunto de datos de imagen; y medir inconsistencias entre el primer y segundo conjuntos de datos de imagen como desplazamientos de fase en el espacio de imagen.

Los aspectos se refieren adicionalmente a la creación de un modelo, del primer y segundo conjuntos de datos de imagen, en donde el modelo puede evaluarse usando la función objetivo.

Los aspectos se refieren adicionalmente, en donde la función objetivo se utiliza para medir un desajuste entre campos de onda de campos de onda sintéticos y el primer y segundo conjuntos de datos de imagen.

Los aspectos se refieren adicionalmente, en donde la función objetivo evalúa un enfoque del modelo.

Los aspectos se refieren adicionalmente, en donde la función objetivo evalúa un grado de semejanza entre el primer y segundo conjuntos de datos de imagen.

Los aspectos de la invención se refieren, por lo tanto, a un sistema que analiza datos de imágenes sísmicas, que comprenden:

un dispositivo de almacenamiento que recibe un primer conjunto de datos de imagen, representando el primer conjunto de datos de imagen datos de imagen obtenidos de receptores sísmicos en una primera posición relativa a una estructura geológica y que recibe un segundo conjunto de datos de imagen, representando el segundo conjunto de datos de imagen datos de imagen obtenidos de receptores sísmicos en una segunda posición relativa a la estructura geológica; y

un módulo de estimación de modelo de velocidad que mide un desplazamiento relativo de los segundos datos de imagen en comparación con los primeros datos de imagen con respecto a la estructura geológica definiendo una función residual, definiendo una función objetivo usando un operador de penalización y calculando un gradiente de la función objetivo.

Los aspectos se refieren adicionalmente, en donde un elemento estructural de la estructura geológica comprende una inmersión estructural.

Los aspectos se refieren adicionalmente, un procesador configurado para comparar el primer conjunto de datos de imagen con el segundo conjunto de datos de imagen y para medir inconsistencias entre el primer y segundo conjuntos de datos de imagen como desplazamientos de fase en el espacio de imagen.

Los aspectos se refieren adicionalmente a, un módulo de función objetivo configurado para construir, del primer y segundo conjuntos de datos, la función objetivo), en donde la función objetivo está configurada para evaluar un modelo creado del primer y segundo conjuntos de datos de imagen.

Los aspectos se refieren adicionalmente, en donde la función objetivo está configurada para medir un desajuste entre campos de onda sintéticos y el primer y segundo conjuntos de datos de imagen.

Los aspectos se refieren adicionalmente, en donde la función objetivo está configurada para evaluar un enfoque del primer y segundo conjuntos de datos de imagen.

Los aspectos se refieren adicionalmente, en donde la función objetivo está configurada para evaluar un grado de semejanza entre el primer y segundo conjuntos de datos de imagen

Estas y otras características y ventajas de esta invención se describen en, o son evidentes a partir de, la siguiente descripción detallada de las realizaciones ilustrativas.

Breve descripción de las figuras

Las realizaciones ilustrativas de la invención se describirán en detalle, con referencia a las siguientes figuras. Debería entenderse que los dibujos no se muestran necesariamente a escala. En ciertos casos, pueden omitirse los detalles que no son necesarios para un entendimiento de la invención o que presentan otros detalles difíciles de percibir. Debería entenderse, por supuesto, que la invención no está necesariamente limitada a las realizaciones particulares ilustradas en el presente documento.

La Figura 1 ilustra un modelo de prueba: una interfaz de densidad horizontal única en un medio de lentitud constante de 0:5 s/km. Los tonos más claros de gris corresponden a valores superiores de densidad.

La Figura 2 ilustra la relación entre el campo de inmersión y el desplazamiento aparente: campos de vector de inmersión (blanco) y de desplazamiento (negro) para diferentes valores de error de lentitud para unos únicos reflectores horizontales:

A m

(a) m =10 %,

A m

(b) m = 0, y

A m

(c ) “ = 10 %.

El valor de lentitud correcto es 0:5 s/km.

La Figura 3 describe el comportamiento de correlaciones locales para el modelo correcto y para uno incorrecto. Cuando el modelo es correcto, el máximo de la correlación de dos imágenes está alineado con la pendiente del reflector (a), de otra manera, el máximo de la correlación se desvía (b).

La Figura 4 ilustra la medida del desplazamiento relativo entre dos imágenes usando correlaciones locales penalizadas. Se consideraron 2 imágenes de experimentos cercanos del reflector en la Figura 2; cada fila (de arriba a abajo) representa un modelo diferente: demasiado alto, correcto y demasiado bajo; las columnas muestran (de izquierda a derecha) el panel de correlación local, el operador de penalización usado y la correlación local penalizada. Para cada modelo, se cogió un punto en el reflector con coordenada lateral x = 3:5 km. La coordenada vertical cambia como una función de parámetros de modelo debido a que la profundidad del reflector cambia.

Obsérvese la asimetría de las correlaciones locales penalizadas para los modelos incorrectos: el valor medio (es decir, la pila a través de los retardos de correlación) de las correlaciones penalizadas nos proporciona una medida del desplazamiento relativo entre dos imágenes.

La Figura 5 ilustra valores de la función objetivo para diferentes modelos en el modelo de lentitud. Se consideraron perturbaciones constantes que variaban del -10 % al 10 % del valor exacto (0,5 s/km).

La Figura 6 ilustra el modelo sencillo que se está considerando en la Figura 1. En este punto, los datos se migran con un modelo de lentitud que es demasiado bajo (a), correcto (b) y demasiado alto (c). El error al cuadrado de la lentitud es constante a través de todo el modelo e igual al 10 % del fondo. El valor de lentitud correcto es 0,5 s/km. La Figura 7 ilustra los tres casos en la Figura 6, donde se calcularon los residuales de acuerdo con la ecuación 2. El signo del residual indica el desplazamiento relativo de una imagen con respecto a la otra.

La Figura 8 ilustra la diferencia entre los residuales de experimentos adyacentes que nos proporciona una estimación de la curvatura de las recopilaciones de la imagen común del dominio de las tomas. Esta información se correlaciona bien con el error de velocidad en el modelo. La Figura 9 ilustra el gradiente de la función objetivo de correlación para un modelo de lentitud que es demasiado bajo (a), correcto (b) y demasiado alto (c). El error al cuadrado de la lentitud es constante en cualquier parte y es igual al 10 % del fondo. La distancia entre las dos tomas es de 60 m. En este punto, se consideraron tres tomas y el gradiente se calculó usando los campos de onda del experimento central.

La Figura 10 ilustra el modelo de lentitud usado para generar los datos acústicos completos. Se han aplicado condiciones de límite absorbente en todos los lados del modelo. Los tonos de gris más claros indican lentitud superior.

La Figura 11 ilustra un modelo de lentitud inicial (a), imagen migrada asociada (b) y recopilaciones de imágenes comunes del dominio de las tomas (c). Obsérvese la curvatura de las recopilaciones, diferentes experimentos representan el reflector a diferentes profundidades.

La Figura 12 ilustra un modelo de lentitud reconstruido (a), imagen migrada final (b) y recopilaciones de imágenes comunes del dominio de las tomas (c). La inversión implica 5 iteraciones de descenso más pronunciado. Obsérvese las recopilaciones horizontales que indican la invarianza de la imagen con respecto a la posición de toma. El área del modelo que está bien iluminada es menor que la extensión real de la perturbación y el modelo invertido tiene un valor superior de lentitud que el medio correcto.

La Figura 13 muestra que la función objetivo se reduce como una función de la iteración y se aplana cuando las recopilaciones de imágenes comunes son horizontales.

La Figura 14 muestra el modelo de lentitud usado para generar los datos acústicos completos. Se han aplicado condiciones de límite absorbente en todos los lados del modelo. Los tonos de gris más claros indican capas más rápidas. La capa superior tiene lentitud de agua (0,66 s/km).

La Figura 15 muestra el modelo de lentitud inicial (a), imagen migrada asociada (b) y recopilaciones de imágenes comunes del dominio de las tomas. Los reflectores están desenfocados y situados incorrectamente de manera grave. Obsérvese la curvatura de las recopilaciones de imágenes comunes que muestran cómo diferentes experimentos representan los reflectores en diferentes profundidades.

La Figura 16 muestra el modelo de lentitud después de 15 interacciones de inversión de forma de onda (a), imagen migrada (b) y recopilaciones de imágenes comunes del dominio de las tomas (c). Obsérvese el enfoque y situación mejorados de los reflectores. Las recopilaciones de imágenes comunes muestran que las imágenes son invariantes en la dirección de toma y el modelo, por lo tanto, es cinemáticamente correcto. Obsérvese el perfil de las diversas capas en el modelo reconstruido.

La Figura 17 muestra la evolución de la función objetivo con iteraciones. Obsérvese la reducción de agudeza de la función objetivo y la rápida convergencia a un modelo cinemáticamente preciso.

La Figura 18 muestra el modelo de lentitud inicial (a), imágenes migradas iniciales (b) y recopilaciones de imágenes comunes del dominio de las tomas (c). El modelo inicial se obtiene suavizando fuertemente y escalando en 0:9 el modelo correcto.

La Figura 19 muestra el modelo actualizado (a), imagen migrada (b) y recopilaciones de imágenes comunes del dominio de las tomas (c) después de 14 iteraciones de descenso más pronunciado.

La Figura 20 muestra la evolución de la función objetivo frente al número de iteraciones.

La Figura 21 ilustra un sistema de tomografía de campo de onda ilustrativo de acuerdo con una realización ilustrativa de esta invención.

La Figura 22 es un diagrama de flujo que ilustra una estimación ilustrativa de un modelo de velocidad de acuerdo con una realización ilustrativa de esta invención.

Descripción detallada

Las realizaciones ilustrativas de esta invención se describirán en relación con el procesamiento y la interpretación de datos, y en particular datos sísmicos. Los sistemas y métodos ilustrativos de esta invención se describirán también en relación con la interpretación y manipulación de datos sísmicos. Sin embargo, para evitar oscurecer innecesariamente la presente invención, la siguiente descripción omite estructuras y dispositivos bien conocidos que pueden mostrarse en forma de diagramas de bloques o resumidos de otra manera.

Para propósitos de explicación, se exponen numerosos detalles para proporcionar un entendimiento minucioso de la presente invención. Sin embargo, debería apreciarse que la presente invención puede ponerse en práctica en una diversidad de maneras más allá de los detalles específicos expuestos en el presente documento.

Adicionalmente, aunque las realizaciones ilustrativas ilustradas en el presente documento muestran los diversos componentes del sistema ubicados conjuntamente, ha de apreciarse que los diversos componentes del sistema pueden ubicarse en porciones distantes de una red distribuida, tal como una red de comunicaciones y/o Internet, o dentro de un sistema seguro, no seguro y/o encriptado especializado. Por lo tanto, debería apreciarse que los componentes del sistema pueden combinarse en uno o más dispositivos o ubicarse conjuntamente en un nodo particular de una red distribuida, tal como una red de comunicaciones. Como se apreciará a partir de la siguiente descripción, y por razones de eficacia computacional, los componentes del sistema pueden disponerse en cualquier ubicación dentro de una red distribuida sin afectar la operación del sistema.

Adicionalmente, debería apreciarse que pueden usarse diversos enlaces para conectar los elementos y pueden ser enlaces alámbricos o inalámbricos, o cualquier combinación de los mismos, o cualquier otro elemento o elementos conocidos o desarrollados más adelante, que puedan suministrar y/o comunicar datos a y desde los elementos conectados. El término módulo como se usa en el presente documento puede hacer referencia a cualquier hardware, software, firmware o combinación de los mismos conocido o desarrollado más adelante que pueda realizar la funcionalidad asociada con ese elemento. Los términos determinar, calcular y computar y variaciones de los mismos, como se usa en el presente documento, se usan de manera intercambiable e incluyen cualquier tipo de metodología, proceso, operación o técnica matemática, que incluye aquellas realizadas por un sistema, tal como un procesador, un sistema experto o red neuronal.

Como se analiza, una realización ilustrativa se refiere al desarrollo de una función objetivo para análisis de velocidad de migración en el dominio de la toma-imagen. La técnica ilustrativa está basada en correlaciones de imagen locales y una reformulación del principio de semejanza que implica únicamente grupos pequeños de imágenes migradas. El criterio de optimización es la minimización del desplazamiento aparente entre las imágenes migradas y conduce de manera natural a un problema de optimización de semejanza diferencial, que se resuelve de manera iterativa usando un método basado en gradiente. La linealidad de los operadores en la función objetivo hace el cálculo del gradiente práctico usando el método de estado adjunto. El enfoque ilustrativo es de onda completa puesto que no se basa en la linealización del procedimiento de extrapolación de onda con respecto a los parámetros de modelo para construir un operador de análisis de velocidad de migración. El sistema resuelve de manera iterativa un problema no lineal con un algoritmo basado en gradiente y la simulación de los campos de onda es completamente no lineal con respecto a los parámetros de modelo. Como se muestra, se prueba el método usando unos pocos ejemplos sintéticos de complejidad creciente.

La Figura 21 ilustra un sistema de tomografía de campo de onda ilustrativo 2100. Además de los componentes bien conocidos, el sistema 2100 comprende un módulo de estimación de modelo de velocidad 2110, un módulo de función objetivo 2120, un controlador de E/S 2130, un controlador/procesador 2140, una memoria 2150, un módulo de formación de imágenes 2160, la pantalla/salida 2170 y el almacenamiento 2180. El sistema 2100 puede estar conectado opcionalmente mediante uno o más enlaces alámbricos o inalámbricos a una o más redes u otros sistemas de tomografía de campo de onda.

En la operación, se reciben datos sísmicos de una o más fuentes y se almacenan en el almacenamiento 2180. A continuación, y en cooperación con el módulo de estimación de modelo de velocidad 2110 y el módulo de función objetivo 2120 (que cooperan con el procesador 2140 y la memoria 2150) se estima el modelo de velocidad utilizando una función objetivo que mide un desajuste entre campos de onda sintéticos y registrados. El gradiente se usa a continuación para actualizar el modelo de velocidad que puede usarse a continuación por el módulo de formación de imágenes 2160 y la pantalla/salida 2170 para emitir imágenes interpretadas.

Más específicamente, las características estructurales en la imagen migrada y la variación de estas características como una función del experimento suministran información valiosa acerca de los parámetros de modelo. En este punto, se muestra que es posible, sin construir recopilaciones de imágenes comunes, extraer información acerca del modelo de un número muy limitado de experimentos usando una relación de deformación entre imágenes migradas. A pesar del hecho de que el sistema puede despreciar la variación del campo de inmersión de la imagen con respecto a los parámetros de modelo, sin embargo, el sistema puede corregir anomalías en el modelo.

Cálculo de la fuente adjunta de la fuente y el receptor

La función objetivo se define como

donde:

es la correlación local de dos imágenes de las tomas cercanas i y (i 1), y P (x, A) es un operador que aniquila el panel de correlación local a lo largo de la dirección de la estructura en cada punto en la imagen.

El método de estado adjunto calcula el gradiente de la función objetivo con respecto a los parámetros de modelo resolviendo dos problemas directos para cada variable de estado del sistema. Este procedimiento resulta ser extremadamente eficaz a nivel computacional puesto que evita perturbar cada parámetro de modelo y calcular la perturbación de variable de estado resultante (Lions, 1972; Plessix, 2006). La solución del problema directo auxiliar requiere diferentes fuentes (fuentes adjuntas) que se calculan diferenciando la función objetivo con respecto a las variables de estado.

Se deriva ahora la expresión matemática para las fuentes adjuntas de la fuente y receptor. El cálculo de la derivada de una funcional con respecto a una función requiere la definición de un producto interno en el espacio funcional apropiado; además, requiere el concepto del diferencial de Fréchet (Vogel, 2002).

Indiquemos las variables de estado del problema, es decir, los campos de onda de fuente y receptor para la toma de orden i, mediante us,i y ur,i, respectivamente. Derivamos la expresión de la fuente adjunta para el campo de onda de fuente. La derivación y la expresión de la fuente adjunta para el campo de onda del receptor son análogas. En primer lugar, hacemos explícita la dependencia del residual n (x) (y por lo tanto de la función objetivo 3 ) en la variable de estado:

donde la integral es una integral triple a través del elemento de volumen dV = dtd^dÁ.

En segundo lugar, calculamos el diferencial de Fréchet de la función objetivo -3 (m ) Con respecto a una perturbación arbitraria de la variable de estado 5us,i. La función objetivo en la ecuación A-1 es la suma a través de las tomas de la energía de los residuales en la ecuación A-2 a través del dominio espacial. Para simplificar la derivación, escribimos en primer lugar el diferencial de la función objetivo como una función del diferencial del residual usando teoría de perturbación:

Expandiendo los cuadrados en la ecuación A-3 y considerando únicamente los términos que son lineales en la perturbación de campo de onda, obtenemos una expresión para el diferencial de Fréchet de la función objetivo:

Los corchetes cuadrados en la ecuación A-4 indican la doble dependencia del diferencial de Fréchet en la variable de estado y su perturbación arbitraria.

Suponiendo que n (x) es diferenciable por Fréchet, calculamos el diferencial de Fréchet usando el límite de la relación incremental:

donde h es un escalar, y 5us,i es una perturbación arbitraria de la variable de estado us,i y pertenece al mismo espacio funcional. El residual n es una funcional lineal de usj, y calculando el límite en la ecuación A-5 obtenemos

Deseamos expresar la ecuación A-6 como un producto interno entre un operador (la derivada de Fréchet) y la perturbación de la variable de estado 5us,¡, realizamos un cambio de variables para hacer que la perturbación de la variable de estado 6us,¡ dependa únicamente de una única variable independiente. Definimos dos nuevas variables „ X

r¡ = ¿;------2 y X = X- el jacobiano de esta transformación es trivial e igual a 1; a continuación, en el nuevo sistema de coordenadas, tenemos

Sri = ¡ P ( x , z ) w ^ x - r i - l p jS u si(Ti,t)-u!, i (Ti,y)RJ(Ti x)dV ', (A -7 )

donde la integral triple es a través del nuevo elemento de volumen dV’ = dtdqdx. A partir de la ecuación A-7, hallamos la expresión del gradiente del residual r, con respecto a la variable de estado us,i.

Combinando la ecuación A-8 en A-4, podemos escribir finamente la ecuación para el gradiente de la función objetivo 3 con respecto a la variable de estado us,i como

donde

son la convolución local y la correlación local del operador de penalización con la imagen a través de la ventana local w, respectivamente. La derivación para la fuente adjunta del lado del receptor es totalmente análoga y conduce a la fórmula

Las ecuaciones A-9 y A-10 son las expresiones para las fuentes adjuntas en las Ecuaciones 10 y 11.

La Figura 22 ilustra un método ilustrativo para tomografía de campo de onda de acuerdo con una realización ilustrativa. Específicamente, el control comienza en la etapa S2200 y continúa a la etapa S2210. En la etapa S2110, se reciben datos sísmicos. A continuación, en la etapa S2120, se estima el modelo de velocidad, que puede usarse más tarde para formación de imágenes de alta calidad de un subsuelo, donde se formula la estimación de velocidad como una función objetivo en la etapa S2130 donde la función objetivo mide un desajuste entre campos de onda sintéticos y registrados. El gradiente determinado se usa a continuación para actualizar el modelo de estimación de velocidad. A continuación, en las etapas S2140 y S2150, se graba el modelo de velocidad actualizado y puede usarse para formación de imágenes de subsuelo. El control a continuación continúa a la etapa S2160 donde finaliza la secuencia de control.

Aunque se ha analizado el diagrama de flujo anteriormente descrito en relación con una secuencia particular de eventos, debería apreciarse que pueden tener lugar cambios en esta secuencia sin afectar de manera material a la operación de la invención. Adicionalmente, la secuencia exacta de eventos no es necesario que tenga lugar como se expone en las realizaciones ilustrativas. Adicionalmente, las técnicas ilustrativas ilustradas en el presente documento no están limitadas a las realizaciones específicamente ilustradas sino que pueden utilizarse también con las otras realizaciones ilustrativas y cada característica descrita puede reivindicarse de manera individual y separada.

Los sistemas, métodos y técnicas de esta invención pueden implementarse en un ordenador de propósito especial, un microprocesador o microcontrolador programado y elemento o elementos de circuitos integrados periféricos, un ASIC u otro circuito integrado, un procesador de señales digitales, un circuito electrónico o lógico de cableado permanente, tal como un circuito de elemento discreto, un dispositivo de lógica programable, tal como PLD, PLA, FPGA, PAL, cualquier medio o similares. En general, cualquier dispositivo que pueda implementar una máquina de estado que a su vez puede implementar la metodología ilustrada en el presente documento puede usarse para implementar los diversos métodos y técnicas de acuerdo con esta invención.

Adicionalmente, los métodos desvelados pueden implementarse fácilmente en software ejecutable por procesador usando entornos de desarrollo de software de objetos u orientados a objetos que proporcionan código fuente portátil que puede usarse en una diversidad de plataformas informáticas o de estación de trabajo. Como alternativa, el sistema desvelado puede implementarse parcial o completamente en hardware usando circuitos lógicos convencionales o diseño de VLSI. Si se usa software o hardware para implementar los sistemas de acuerdo con esta invención es dependiente de los requisitos de velocidad y/o eficiencia del sistema, la función particular y los sistemas de software o hardware particulares o los sistemas de microprocesadores o microordenadores informáticos que se estén utilizando. Los sistemas, métodos y técnicas ilustrados en el presente documento pueden implementarse fácilmente en hardware y/o software usando cualquier sistema o estructuras, dispositivos y/o software conocidos o desarrollados más adelante por los expertos en la materia en la técnica aplicable de la descripción funcional proporcionada en el presente documento y con un conocimiento básico general de las artes de la informática y la geología.

Además, los métodos desvelados pueden implementarse fácilmente en software que puede almacenarse en un medio de almacenamiento legible por ordenador no transitorio, ejecutarse en un ordenador de propósito general programado con la cooperación de un controlador y memoria, un ordenador de propósito especial, un microprocesador o similares. Los sistemas y métodos de esta invención pueden implementarse como un programa embebido en ordenador personal tal como una miniaplicación, JAVA® o CGI script, en C o C++, Fortran, o similares, como un recurso que reside en un servidor o estación de trabajo informática, como una rutina embebida en un sistema especializado o componente de sistema o similares. El sistema puede implementarse también físicamente incorporando el sistema y/o método en un sistema de software y/o hardware, tal como los sistemas de hardware y software de un dispositivo de interpretación sísmica especializado.

Referencias:

Un experto en la materia tendría conocimiento de las siguientes referencias que se incorporan en el presente documento por referencia en su totalidad:

Al-Yahya, K., 1989, Velocity analysis by iterative profile migration: Geophysics, 54, págs. 718-729.

Albertin, U., P. Sava, J. Etgen y M. Maharramov, 2006, Adjoint wave-equation velocity analysis: presentado en 74th Ann. Internat. Mtg., SEG.

Baysal, E., D. D. Kosloff y J. W. C. Sherwood, 1983, Reverse time migration: Geophysics, 48, 1514-1524. Biondi, B. y P. Sava, 1999, Wave-equation migration velocity analysis: presentado en 69th Ann. Internat. Mtg., SEG. Biondi, B. y W.W. Symes, 2004, Angle-domain commonimage gathers for migration velocity analysis by wavefield continuation imaging: Geophysics, 69, 1283-1298.

Bishop, T. N., K. P. Bube, R. T. Cutler, R. T. Langan, P. L. Love, J. R. Resnick, R. T. Shuey, D. A. Spindler y H. W. Wyld, 1985, Tomographic determination of velocity and depth in laterally varying media: Geophysics, 50, 903-923. Carcione, J., 2007, Wave fields in real media wave propagation in anisotropic, anelastic, porous y electromagnetic media: Elsevier. Handbook of Geophysical Exploration.

Chang, W., P. Carrion y G. Beylkin, 1987, Wavefront sets of solutions to linearised inverse scattering problems: Inverse Problems, 3, 683 - 690.

Chavent, G. y C. A. Jacewitz, 1995, Determination of background velocities by multiple migration fitting: Geophysics, 60, 476-490.

Claerbout, J. F., 1985, Imaging the earth's interior: Blackwell Publishing.

Claerbout, J. F. y F. Muir, 1973, Robust modeling with erratic data: Geophysics, 38, 826-844.

de Vries, D. y A. J. Berkhout, 1984, Velocity analysis based on minimum entropy: Geophysics, 49, 2132-2142. Farmer, P. A., I. F. Jones, H. Zhou, R. I. Bloor y M. C. Goodwin, 2006, Application of reverse time migration to complex imaging problems: First Break, 24, 65-73.

Faye, J.-P. y J.-P. Jeannot, 1986, Prestack migration velocity analysis from focusing depth imaging: presentado en 56th Ann. Internat. Mtg., Soc. of Expl. Geophys.

Fomel, S., 2002, Applications of plane-wave destruction filters: Geophysics, 67, 1946-1960.

Fowler, P., 1985, Migration velocity analysis by optimization: Technical Report 44, Stanford Exploration Project. Gazdag, J. y P. Sguazzero, 1984, Migration of seismic data by phase shift plus interpolation: Geophysics, 49. Hale, D., 2006, An efficient method for computing local cross-correlations of multidimensional signals: Technical Report CWP-544, Center for Wave Phenomena, Colorado School of Mines.

-, 2007a, Local dip filtering with directional laplacian: Technical Report CWP-567, Center for Wave Phenomena, Colorado School of Mines.

-, 2007b, A method for estimating apparent displacement vectors from time-lapse seismic data: T echnical Report CWP-566, Center for Wave Phenomena, Colorado School of Mines.

Kelly, S., J. R. Martinez, B. tsimelzon y S. Crawley, 2010, A comparison of inversion results for two full-waveform methods that utilize the lowest frequencies in dual-sensorrecordings: presentado en 80th Ann. Internat. Mtg., Soc. of Expl. Geophys.

Lambar'e, G., M. Alerini, R. Baina y P. Podvin, 2004, Stereotomography: a semiautomatic approach for velocity macromodel estimation: Geophysical Prospecting, 52, 671-681.

Lions, J., 1972, Nonhomogeneous boundary value problems and applications: Springer Verlag, Berlín.

Liu, F., D. W. Hanson, N. D. Whitmore, R. S. Day y R. H. Stolt, 2006, Toward a unified analysis for source plane wave migration: Geophysics, 71, S129-S139.

Ma, Y., D. Hale, Z. Meng y B. Gong, 2010, Full waveform inversion with image-guided gradient: presentado en 80th Ann. Internat. Mtg., Soc. of Expl. Geophys. McMechan, G. A., 1983, Migration by extrapolación of time-dependent boundary values: Geophysical Prospecting, 31,413-420.

Morton, S. A. y C. C. Ober, 1998, Faster shot-record migration using phase encoding: presentado en 68th Ann. Internat. Mtg., SEG.

Mulder, W. A. y A. P. E. Kroode, 2002, Automatic velocity analysis by differential semblance optimization: Geophysics, 67, 1184 - 1191.

Perrone, F. y P. Sava, 2011,Wave-equation migration with dithered plane waves: Geophysical Prospecting. Perrone, F., P. Sava, C. Andreoletti y N. Bienati, 2012, Wave-equation migration velocity analysis by image warping: Technical Report CWP-708, Center for Wave Phenomena, Colorado School of Mines.

Plessix, R.-E., 2006, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications: Geophys. J. Int., 167, 495-503.

Pratt, R. G., 1999, Seismic waveform inversion in the frequency domain, part 1: Theory and verification in a physical scale model: Geophysics, 64, 888-901.

Rickett, J. y P. Sava, 2002, Offset and angle-domain common image-point gathers for shot-profile migration: Geophysics, 67, 883-889.

Romero, L. A., D. C. Ghiglia, C. C. Ober y S. A. Morton, 2000, Phase encoding of shot records in prestack migration: Geophysics, 65, 426-436.

Santosa, F. y W. W. Symes, 1989, An analysis of leastsquares velocity inversion, in Geophysical Monograph n° 4: Soc. of Expl. Geophys.

Sava, P. y B. Biondi, 2004, Wave-equation migration velocity analysis, i. theory: Geophysical Prospecting, 52, 593 606.

Sava, P., B. Biondi y J. Etgen, 2005, Wave-equation migration velocity analysis by focusing diffractions and reflections:Geophysics, 70, U19-U27.

Sava, P. y S. Fomel, 2003, Angle-domain common-image gathers by wavefield continuation methods: Geophysics, 68, 1065-1074.

Sava, P. e I. Vasconcelos, 2009, Extended commonimage-point gathers for wave-equation migration: presentado en 71th Ann. Internat. Mtg., European Association of Geoscientists and Engineers.

Shen, P. y W. W. Symes, 2008, Automatic velocity analysis via shot profile migration: Geophysics, 73, VE49-VE59. Sirgue, L., O.I. Barkved, J. Dellinger, J. Etgen, U. Albertin y J. H. Kommedal, 2010, Full waveform inversion: the next leap forward in imaging at valhall: First Break, 28, 65-70.

Sirgue, L. y R. Pratt, 2004, Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies: Geophysics, 69, 231-248.

Soubaras, R., 2006, Modulated-shot migration: 76th Ann. Internat. Mtg., Soc. of Expl. Geophys., 2430-2433. Stoffa, P. L., J. T. Fokkema, R. M. de Luna Freire y W. P. Kessiger, 1990, Split-step fourier migration: Geophysics, 55. Stoffa, P. L., M. K. Sen, R. K. Seifoullaev, R. C. Pestana y J. T. Fokkema, 2006, Plane-wave depth migration: Geophysics, 71, S261-S272.

Stolt, R. H., 1978, Migration by fourier transform: Geophysics, 43.

Symes, W. W., 1991, Layered velocity inversion: A model problem from reflection seismology problem from reflection seismology: SIAM J. on Math. Analysis, 22, págs. 680-716. -, 2008, Migration velocity analysis and waveform inversion: Geophysical Prospecting, 56, 765-790.

Symes, W. W. y J. J. Carazzone, 1991, Velocity inversion by differential semblance optimization: Geophysics, 56, 654-663.

Taner, M. T. y F. Koehler, 1969, Velocity spectra - digital computer derivation and applications of velocity functions: Geophysics, 34, 859 - 881.

Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49, 71-92. -, 2005, Inverse problem theory and methods for modelparameter estimation: SIAM.

van Leeuwen, T. y W. A. Mulder, 2008, Velocity analysis based on data correlation: Geophysical Prospecting, 56, 791- 803.

van Vliet, L. J. y P. W. Verbeek, 1995, Estimators for orientation and anisotropy in digitalized images: Proceedings of the first annual conference of the Advanced School for Computing and Imaging, ASCI'95, 442-450.

Vogel, C., 2002, Computational methods for inverse problems: SIAM.

Whitmore, N., 1995, An imaging hierarchy for commonangle plane wave seismograms: PhD thesis, University of Tulsa.

Woodward, M. J., 1992, Wave-equation tomography: Geophysics, 57, 15-26.

Xie, X. B. y H. Yang, 2008, The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis: Geophysics, 73, S241-S249.

Yang, T. y P. Sava, 2010, Wave-equation migration velocity analysis with extended common-image-point gathers: presentado en 80th Ann. Internat. Mtg., Soc. of Expl. Geophys. -, 2011,Wave-equation migration velocity analysis with time-shift imaging: Geophysical Prospecting, 59, 635-650.

Yilmaz, O. y R. Chambers, 1984, Migration velocity analysis by wave-field extrapolation: Geophysics, 49, P. 1664 1674.

Zhang, Y., J. Sun, C.Notfors, S. H. Gray, L. Chernis y J. Young, 2005, Delayed-shot 3d depth migration: Geophysics, 70, E21-E28.

Claims

REIVINDICACIONES

1. Un método implementado por ordenador (S2100) de análisis de datos de imágenes sísmicas, comprendiendo el método:

recibir un primer conjunto de datos de imagen (S2110), representando el primer conjunto de datos de imagen datos de imagen de una estructura geológica obtenidos de receptores sísmicos en una primera posición relativa a la estructura geológica;

recibir un segundo conjunto de datos de imagen (S2110), representando el segundo conjunto de datos de imagen datos de imagen de la estructura geológica obtenidos de receptores sísmicos en una segunda posición relativa a la estructura geológica; y caracterizado por que el método comprende adicionalmente

medir un desplazamiento relativo (S2130) de los segundos datos de imagen en comparación con los primeros datos de imagen con respecto a la estructura geológica comparando directamente el primer conjunto de datos de imagen con el segundo conjunto de datos de imagen definiendo una función residual, definiendo una función objetivo usando un operador de penalización y calculando un gradiente de la función objetivo.

2. El método de la reivindicación 1, en donde un elemento estructural de la estructura geológica comprende una inmersión estructural.

3. El método de la reivindicación 1, en donde medir un desplazamiento relativo comprende:

comparar, usando un ordenador (2100), el primer conjunto de datos de imagen con el segundo conjunto de datos de imagen; y

medir inconsistencias (S2130) entre el primer y segundo conjuntos de datos de imagen como desplazamientos de fase en el espacio de imagen.

4. El método de la reivindicación 3, en donde medir un desplazamiento relativo comprende adicionalmente crear un modelo, del primer y segundo conjuntos de datos de imagen (S2120), en donde el modelo puede evaluarse usando la función objetivo (S2130).

5. El método de la reivindicación 4, en donde la función objetivo (S2130) se utiliza para medir un desajuste entre campos de onda sintéticos y el primer y segundo conjuntos de datos de imagen.

6. El método de la reivindicación 4, en donde la función objetivo evalúa un enfoque del modelo.

7. El método de la reivindicación 4, en donde la función objetivo evalúa un grado de semejanza entre el primer y segundo conjuntos de datos de imagen.

8. Un sistema (S2100) que analiza datos de imágenes sísmicas, que comprende:

un dispositivo de almacenamiento (2180) que comprende un primer conjunto de datos de imagen, representando el primer conjunto de datos de imagen datos de imagen (S2110) obtenidos de receptores sísmicos en una primera posición relativa a una estructura geológica y que recibe un segundo conjunto de datos de imagen (S2110), representando el segundo conjunto de datos datos de imagen obtenidos de receptores sísmicos en una segunda posición relativa a la estructura geológica; y caracterizado por comprender adicionalmente

un módulo de estimación de modelo de velocidad (2110) configurado para medir un desplazamiento relativo (S2130) de los segundos datos de imagen en comparación con los primeros datos de imagen con respecto a la estructura geológica definiendo una función residual, definiendo una función objetivo usando un operador de penalización y calculando un gradiente de la función objetivo.

9. El sistema de la reivindicación 8, en donde un elemento estructural de la estructura geológica comprende una inmersión estructural.

10. El sistema de la reivindicación 8, que comprende adicionalmente:

un procesador (2140) configurado para comparar el primer conjunto de datos de imagen con el segundo conjunto de datos de imagen y para medir inconsistencias entre el primer y segundo conjuntos de datos de imagen como desplazamientos de fase en el espacio de imagen.

11. El sistema de la reivindicación 10, que comprende adicionalmente un módulo de función objetivo (2120) configurado para construir, del primer y segundo conjuntos de datos de imagen, la función objetivo (S2130), en donde la función objetivo está configurada para evaluar un modelo creado del primer y segundo conjuntos de datos de imagen.

12. El sistema de la reivindicación 11, en donde la función objetivo está configurada para medir un desajuste (S2130) entre campos de onda sintéticos y el primer y segundo conjuntos de datos de imagen.

13. El sistema de la reivindicación 11, en donde la función objetivo está configurada para evaluar un enfoque del primer y segundo conjuntos de datos de imagen.

14. El sistema de la reivindicación 11, en donde la función objetivo está configurada para evaluar un grado de semejanza entre el primer y segundo conjuntos de datos de imagen.

15. Un medio de almacenamiento de información legible por ordenador no transitorio (2150) que tiene almacenado en el mismo instrucciones, que cuando se ejecutan por uno o más ordenadores (2100), realizan las etapas de la reivindicación 1.