Inversión sísmica lineal

La pieza central de la inversión lineal sísmica es el modelo directo que modela la generación de los datos experimentales recopilados. ^{[ 1 ]} Según Wiggins (1972), ^{[ 4 ]} proporciona una relación funcional (computacional) entre los parámetros del modelo y los valores calculados para las trazas observadas. Dependiendo de los datos sísmicos recopilados, este modelo puede variar desde las ecuaciones de onda clásicas para predecir el desplazamiento de partículas o la presión del fluido para la propagación de ondas sonoras a través de rocas o fluidos, hasta algunas variantes de estas ecuaciones clásicas. Por ejemplo, el modelo directo en Tarantola (1984) ^{[ 5 ]} es la ecuación de onda para la variación de presión en un medio líquido durante la propagación de ondas sísmicas , mientras que, al asumir capas de velocidad constante con interfaces planas, Kanasewich y Chiu (1985) ^{[ 6 ]} utilizaron el modelo de braquistótrona de John Bernoulli para el tiempo de viaje de un rayo a lo largo de una trayectoria. En Cooke y Schneider (1983), ^{[ 1 ]} el modelo es un algoritmo de generación de trazas sintéticas expresado como en la ecuación. 3, donde R(t) se genera en el dominio Z mediante una fórmula recursiva. Cualquiera que sea la forma en que aparezca el modelo directo, es importante que no solo prediga los datos de campo recopilados, sino que también modele cómo se generan dichos datos. Por lo tanto, el modelo directo de Cooke y Schneider (1983) ^{[ 1 ]} solo puede utilizarse para invertir datos CMP, ya que el modelo invariablemente asume que no hay pérdida por propagación al imitar la respuesta de una tierra lateralmente homogénea a una fuente de onda plana.

1. ${\displaystyle t=\sum _{i=1}^{n}{\frac {{\big [}(x_{i}-x_{i-1})^{2}+(y_{i}-y_{i-1})^{2}+(z_{i}-z_{i-1})^{2}{\big ]}^{\frac {1}{2}}}{v_{i}}}}$
donde t es el tiempo de viaje del rayo, x, y, z son las coordenadas de profundidad y vi es la velocidad constante entre las interfaces i − 1 e i.
2. ${\displaystyle \left[{\frac {1}{K({\vec {r}})}}{\frac {\partial ^{2}}{\partial t^{2}}}-\nabla \cdot {\big (}{\frac {1}{\rho ({\vec {r}}{\big )}}}\nabla )\right]U({\vec {r}},t)=s({\vec {r}},t)}$
dónde${\displaystyle K({\vec {r}})}$representan el módulo de compresibilidad,${\displaystyle \rho ({\vec {r}})}$densidad,${\displaystyle s({\vec {r}},t)}$la fuente de las ondas acústicas y${\displaystyle U({\vec {r}},t)}$ la variación de presión.
3. ${\displaystyle s(t)=w(t)*R(t)}$

donde s ( t ) = traza sintética, w ( t ) = ondícula fuente y R ( t ) = función de reflectividad.

Un proceso numérico importante en la modelización inversa consiste en minimizar la función objetivo, definida como la diferencia entre los datos sísmicos de campo recopilados y los datos sísmicos calculados numéricamente. Las funciones objetivo clásicas incluyen la suma de las desviaciones cuadráticas entre los datos experimentales y numéricos, como en los métodos de mínimos cuadrados , la suma de la magnitud de la diferencia entre los datos de campo y numéricos, o alguna variante de estas definiciones. Independientemente de la definición utilizada, la solución numérica del problema inverso se obtiene como un modelo terrestre que minimiza la función objetivo.

Además de la función objetivo, en el procedimiento de modelado inverso también se incorporan otras restricciones, como parámetros de modelo conocidos e interfaces de capas conocidas en algunas regiones de la Tierra. Estas restricciones, según Francis (2006) ^{[ 3 ],} ayudan a reducir la no unicidad de la solución de inversión al proporcionar información a priori que no está contenida en los datos invertidos, mientras que Cooke y Schneider (1983) ^{[ 1 ]} informan de su utilidad para controlar el ruido y cuando se trabaja en un área geofísicamente bien conocida.

El objetivo del análisis matemático del modelado inverso es convertir el problema inverso lineal generalizado en un álgebra matricial simple considerando todos los componentes descritos en secciones anteriores, a saber: modelo directo, función objetivo, etc. Generalmente, los datos sísmicos generados numéricamente son funciones no lineales de los parámetros del modelo terrestre. Para eliminar la no linealidad y crear una plataforma para la aplicación de conceptos de álgebra lineal , el modelo directo se linealiza mediante una expansión utilizando una serie de Taylor, como se muestra a continuación. Para más detalles, véase Wiggins (1972), ^{[ 4 ]} Cooke y Schneider (1983). ^{[ 1 ]}

Consideremos un conjunto de${\displaystyle m}$observaciones de campo sísmico${\displaystyle F_{j}}$, para${\displaystyle j=1,\ldots ,m}$y un conjunto de${\displaystyle n}$parámetros del modelo terrestre${\displaystyle p_{i}}$ser invertido para, para${\displaystyle i=1,\ldots ,n}$Las observaciones de campo pueden representarse en cualquiera de los siguientes formatos:${\displaystyle {\vec {F}}\,({\vec {p}})}$o${\displaystyle F_{j}\,(p_{i})}$, dónde${\displaystyle {\vec {p}}}$y${\displaystyle {\vec {F}}\,({\vec {p}})}$son representaciones vectoriales de los parámetros del modelo y las observaciones de campo en función de los parámetros terrestres. De manera similar, para${\displaystyle q_{i}}$representando estimaciones de los parámetros del modelo,${\displaystyle {\vec {F}}\,({\vec {q}})}$es el vector de datos sísmicos calculados numéricamente utilizando el modelo directo de la Sec. 1.3. Expansión en serie de Taylor de${\displaystyle {\vec {F}}\,({\vec {p}})}$acerca de${\displaystyle {\vec {q}}}$Se proporciona a continuación.

4. ${\displaystyle {\vec {F}}\,({\vec {p}})={\vec {F}}\,({\vec {q}})+({\vec {p}}-{\vec {q}}){\frac {\partial {\vec {F}}\,({\vec {q}})}{\partial {\vec {p}}}}+({\vec {p}}-{\vec {q}})^{2}{\frac {\partial ^{2}{\vec {F}}\,({\vec {q}})}{\partial {\vec {p}}^{2}}}+O({\vec {p}}-{\vec {q}})^{3}}$
Al linealizar eliminando los términos no lineales (términos con (p⃗ − ⃗q) de orden 2 o superior), la ecuación se convierte en:
5. ${\displaystyle {\vec {F}}\,({\vec {p}})-{\vec {F}}\,({\vec {q}})=({\vec {p}}-{\vec {q}}){\frac {\partial {\vec {F}}\,({\vec {q}})}{\partial {\vec {p}}}}}$
Considerando que${\displaystyle {\vec {F}}}$tiene${\displaystyle m}$componentes y${\displaystyle {\vec {p}}}$y${\displaystyle {\vec {q}}}$tener${\displaystyle n}$componentes, la forma discreta de la ecuación 5 da como resultado un sistema de${\displaystyle m}$ecuaciones lineales en${\displaystyle n}$variables cuya forma matricial se muestra a continuación.
6. ${\displaystyle \Delta {\vec {F}}=\mathbf {A} \,\Delta {\vec {p}}}$
7. ${\displaystyle \Delta {\vec {F}}={\begin{bmatrix}F_{1}({\vec {p}})-F_{1}({\vec {q}})\\\vdots \\F_{m}({\vec {p}})-F_{m}({\vec {q}})\end{bmatrix}}}$
8. ${\displaystyle \Delta {\vec {p}}={\vec {p}}-{\vec {q}}={\begin{bmatrix}p_{1}-q_{1}\\\vdots \\p_{n}-q_{n}\end{bmatrix}}}$
9. ${\displaystyle \mathbf {A} ={\begin{bmatrix}{\frac {\partial F_{1}({\vec {q}})}{\partial p_{1}}}&{\frac {\partial F_{1}({\vec {q}})}{\partial p_{2}}}&\cdots &{\frac {\partial F_{1}({\vec {q}})}{\partial p_{n}}}\\{\frac {\partial F_{2}({\vec {q}})}{\partial p_{1}}}&\cdots &{\frac {\partial F_{2}({\vec {q}})}{\partial p_{n-1}}}&{\frac {\partial F_{2}({\vec {q}})}{\partial p_{n}}}\\\vdots &{\frac {\partial F_{j}({\vec {q}})}{\partial p_{i}}}&\vdots &\vdots \\{\frac {\partial F_{m}({\vec {q}})}{\partial p_{1}}}&{\frac {\partial F_{m}({\vec {q}})}{\partial p_{2}}}&\cdots &{\frac {\partial F_{m}({\vec {q}})}{\partial p_{n}}}\\\end{bmatrix}}}$

${\displaystyle \Delta {\vec {F}}}$se denomina vector de diferencias en Cooke y Schneider (1983). ^{[ 1 ]} Tiene un tamaño de${\displaystyle m\times 1}$y sus componentes son la diferencia entre el rastro observado y los datos sísmicos calculados numéricamente.${\displaystyle \Delta {\vec {p}}}$es el vector corrector de tamaño${\displaystyle n\times 1}$, mientras${\displaystyle \mathbf {A} }$se denomina matriz de sensibilidad. Tiene un tamaño de${\displaystyle m\times n}$y sus comentarios son tales que cada columna es la derivada parcial de un componente de la función directa con respecto a uno de los parámetros desconocidos del modelo terrestre. De manera similar, cada fila es la derivada parcial de un componente de la traza sísmica calculada numéricamente con respecto a todos los parámetros desconocidos del modelo.

${\displaystyle {\vec {F}}\,({\vec {q}})}$se calcula a partir del modelo directo, mientras que${\displaystyle {\vec {F}}\,({\vec {p}})}$son los datos experimentales. Por lo tanto,${\displaystyle \Delta {\vec {F}}}$es una cualidad conocida. Por otro lado,${\displaystyle \Delta {\vec {p}}}$es desconocido y se obtiene mediante la solución de la ecuación 10. Esta ecuación es teóricamente resoluble solo cuando${\displaystyle \mathbf {A} }$es invertible, es decir, si es una matriz cuadrada de modo que el número de observaciones${\displaystyle m}$es igual al número${\displaystyle n}$de parámetros terrestres desconocidos. Si este es el caso, el vector corrector desconocido${\displaystyle \Delta {\vec {p}}}$, se resuelve como se muestra a continuación, utilizando cualquiera de los solucionadores clásicos directos o iterativos para la solución de un conjunto de ecuaciones lineales.

10. ${\displaystyle \Delta {\vec {p}}=\mathbf {A} ^{-1}\,\Delta {\vec {F}}}$

En la mayoría de las aplicaciones de inversión sísmica , hay más observaciones que el número de parámetros terrestres que se deben invertir, es decir${\displaystyle m>n}$, lo que conduce a un sistema de ecuaciones matemáticamente sobredeterminado. Como resultado, la ecuación 10 no es teóricamente resoluble y no se puede obtener una solución exacta. ^{[ 6 ]} Se obtiene una estimación del vector corrector utilizando el procedimiento de mínimos cuadrados para encontrar el vector corrector.${\displaystyle \Delta {\vec {p}}}$que minimiza${\displaystyle {\vec {e}}\,^{T}{\vec {e}}}$, que es la suma de los cuadrados del error,${\displaystyle {\vec {e}}}$. ^{[ 6 ]}

El error${\displaystyle {\vec {e}}}$es dado por

11. ${\displaystyle {\vec {e}}=\Delta {\vec {F}}-\mathbf {A} \,\Delta {\vec {p}}}$

En el procedimiento de mínimos cuadrados, el vector corrector que minimiza${\displaystyle {\vec {e}}\,^{T}{\vec {e}}}$se obtiene como se muestra a continuación.

12. ${\displaystyle {\begin{aligned}\mathbf {A} \,\Delta {\vec {p}}&=\Delta {\vec {F}}\\\mathbf {A} ^{T}\mathbf {A} \,\Delta {\vec {p}}&=\mathbf {A} ^{T}\Delta {\vec {F}}\end{aligned}}}$

De este modo,

13. ${\displaystyle \Delta {\vec {p}}=(\mathbf {A} ^{T}\mathbf {A} )^{-1}\,\mathbf {A} ^{T}\Delta {\vec {F}}}$

De las discusiones anteriores, la función objetivo se define como la${\displaystyle L_{1}}$o${\displaystyle L_{2}}$norma de${\displaystyle \Delta {\vec {p}}}$dado por ${\displaystyle \sum _{j=0}^{n}|\Delta p_{j}|}$o${\displaystyle \sum _{j=0}^{n}|\Delta p_{j}|^{2}}$o de${\displaystyle \Delta {\vec {F}}}$dado por${\displaystyle \sum _{i=0}^{m}|\Delta F_{i}|}$o${\displaystyle \sum _{i=0}^{m}|\Delta F_{i}|^{2}}$.

El procedimiento generalizado para invertir cualquier dato sísmico experimental para${\displaystyle m=n}$o${\displaystyle m>n}$, utilizando la teoría matemática para el modelado inverso, como se describió anteriormente, se muestra en la Fig. 1 y se describe a continuación.

Se proporciona una estimación inicial de la impedancia del modelo para iniciar el proceso de inversión. El modelo directo utiliza esta estimación inicial para calcular datos sísmicos sintéticos, que se restan de los datos sísmicos observados para calcular el vector de diferencia.

1. Una estimación inicial de la impedancia del modelo${\displaystyle {\vec {q}}}$Se proporciona para iniciar el proceso de inversión.
2. Un conjunto de datos sísmicos sintéticos${\displaystyle {\vec {F}}({\vec {q}})}$se calcula mediante el modelo directo, utilizando la impedancia del modelo anterior.
3. El vector de diferencias${\displaystyle {\vec {F}}({\vec {p}})-{\vec {F}}({\vec {q}})}$se calcula como la diferencia entre los datos sísmicos experimentales y sintéticos.
4. La matriz de sensibilidad${\displaystyle \mathbf {A} }$se calcula en este valor del perfil de impedancia.
5. Usando${\displaystyle \mathbf {A} }$y el vector de diferencia del 3 anterior, el vector corrector${\displaystyle \Delta {\vec {p}}}$se calcula. Se obtiene un nuevo perfil de impedancia como
   14. ${\displaystyle {\vec {p}}={\vec {q}}+\Delta {\vec {p}}}$
6. El${\displaystyle L_{1}}$o${\displaystyle L_{2}}$La norma del vector corrector calculado se compara con un valor de tolerancia proporcionado. Si la norma calculada es menor que la tolerancia, el procedimiento numérico concluye y el perfil de impedancia invertido para la región terrestre viene dado por${\displaystyle {\vec {p}}}$de la ecuación 14. Por otro lado, si la norma es mayor que la tolerancia, las iteraciones a través de los pasos 2-6 se repiten pero con un perfil de impedancia actualizado como se calcula a partir de la ecuación 14. La figura 2 ^{[ 7 ]} muestra un ejemplo típico de actualización del perfil de impedancia durante el proceso de iteración sucesiva. Según Cooke y Schneider (1983), ^{[ 1 ]} el uso de la estimación corregida de la ecuación 14 como nueva estimación inicial durante la iteración reduce el error.

Fuente: ^{[ 8 ]}

Comenzamos con un ejemplo para invertir los valores de los parámetros terrestres a partir de la distribución de temperatura en función de la profundidad en una región terrestre determinada. Si bien este ejemplo no se relaciona directamente con la inversión sísmica, ya que no intervienen ondas acústicas viajeras, introduce una aplicación práctica de la técnica de inversión de manera comprensible, antes de pasar a las aplicaciones sísmicas. En este ejemplo, la temperatura de la Tierra se mide en puntos discretos de un pozo mediante la colocación de sensores de temperatura a las profundidades objetivo. Suponiendo un modelo directo de distribución lineal de la temperatura con la profundidad, se invierten dos parámetros a partir de las mediciones de temperatura en función de la profundidad.

El modelo directo viene dado por

15. ${\displaystyle {\vec {F}}({\vec {q}})={\vec {T}}=a+bz}$

dónde${\displaystyle {\vec {q}}=[a,b]}$. Por lo tanto, la dimensión de${\displaystyle {\vec {q}}}$es 2, es decir, el número de parámetros invertidos es 2.

El objetivo de este algoritmo de inversión es encontrar${\displaystyle {\vec {p}}}$, que es el valor de${\displaystyle [a,b]}$que minimiza la diferencia entre la distribución de temperatura observada y las obtenidas utilizando el modelo directo de la ecuación 15. Considerando la dimensión del modelo directo o el número de observaciones de temperatura a ser${\displaystyle n}$, los componentes del modelo directo se escriben como

16. ${\displaystyle {\begin{aligned}T_{1}&=a+bz_{1}\\T_{2}&=a+bz_{2}\\\vdots \\T_{n-1}&=a+bz_{n-1}\\T_{n}&=a+bz_{n}\\\end{aligned}}}$
de modo que${\displaystyle {\vec {F}}({\vec {q}})=T}$
17. ${\displaystyle \mathbf {A} ={\begin{bmatrix}1&z_{1}\\1&z_{2}\\\vdots &\vdots \\1&z_{n-1}\\1&z_{n}\\\end{bmatrix}}}$

Presentamos los resultados de Marescot (2010) ^{[ 8 ]} para el caso de${\displaystyle n=2}$para los cuales los valores de temperatura observados a profundidades fueron${\displaystyle T_{1}=19^{\circ }C}$en${\displaystyle z=2m}$y${\displaystyle T_{2}=22^{\circ }C}$en${\displaystyle z=8m}$Estos datos experimentales se invirtieron para obtener valores de parámetros terrestres de${\displaystyle a=0.5}$y${\displaystyle b=18^{\circ }C}$. Para un caso más general con un gran número de observaciones de temperatura, la Fig. 4 muestra el modelo lineal directo final obtenido al usar los valores invertidos de${\displaystyle a}$y${\displaystyle b}$La figura muestra una buena concordancia entre los datos experimentales y numéricos.

Fuente: ^{[ 8 ]}

Este ejemplo invierte la velocidad de la capa terrestre a partir de los tiempos de llegada de las ondas sísmicas registradas. La figura 5 muestra las estimaciones iniciales de velocidad y los tiempos de llegada registrados en campo, mientras que la figura 6a muestra el modelo de velocidad heterogéneo invertido , que es la solución del algoritmo de inversión obtenida tras 30 iteraciones . Como se observa en la figura 6b, existe una buena concordancia entre los tiempos de llegada finales obtenidos a partir del modelo directo utilizando la velocidad invertida y los tiempos de llegada registrados en campo. Utilizando estas soluciones, se reconstruyó la trayectoria del rayo, la cual se muestra como altamente tortuosa a través del modelo terrestre, tal como se observa en la figura 7.

Este ejemplo, tomado de Cooke y Schneider (1983), ^{[ 1 ]} muestra la inversión de una traza sísmica CMP para obtener el perfil de impedancia del modelo terrestre (producto de la densidad y la velocidad). La traza sísmica invertida se muestra en la Fig. 8, mientras que la Fig. 9a muestra el perfil de impedancia invertido con la impedancia inicial de entrada utilizada para el algoritmo de inversión. Junto con la traza sísmica, también se registró un registro de impedancia de la región terrestre, como se muestra en la Fig. 9b. Las figuras muestran una buena concordancia entre el registro de impedancia registrado y la impedancia numérica invertida a partir de la traza sísmica.