Problema inverso

La primera dificultad refleja un problema crucial: nuestras observaciones no contienen suficiente información y se requieren datos adicionales. Estos datos adicionales pueden provenir de información física previa sobre los valores de los parámetros, su distribución espacial o, de forma más general, su interdependencia. También pueden provenir de otros experimentos: por ejemplo, podemos considerar la integración de datos registrados por gravímetros y sismógrafos para una mejor estimación de las densidades. La integración de esta información adicional es, fundamentalmente, un problema de estadística . Esta disciplina es la que puede responder a la pregunta: ¿Cómo combinar cantidades de distinta naturaleza? Seremos más precisos en la sección «Enfoque bayesiano» que se presenta a continuación.

En lo que respecta a los parámetros distribuidos, la información previa sobre su distribución espacial a menudo consiste en información sobre algunas derivadas de estos parámetros distribuidos. Además, es práctica común, aunque algo artificial, buscar el modelo "más simple" que se ajuste razonablemente a los datos. Esto generalmente se logra penalizando la${\displaystyle L^{1}}$norma del gradiente (o la variación total ) de los parámetros (este enfoque también se conoce como maximización de la entropía). También se puede simplificar el modelo mediante una parametrización que introduzca grados de libertad solo cuando sea necesario.

También se puede integrar información adicional mediante restricciones de desigualdad en los parámetros del modelo o algunas funciones de ellos. Dichas restricciones son importantes para evitar valores poco realistas para los parámetros (valores negativos, por ejemplo). En este caso, el espacio generado por los parámetros del modelo ya no será un espacio vectorial, sino un subconjunto de modelos admisibles denotado por${\displaystyle P_{\text{adm}}}$en la secuela.

Como se mencionó anteriormente, el ruido puede ser tal que nuestras mediciones no sean la imagen de ningún modelo, por lo que no podemos buscar un modelo que produzca los datos, sino más bien buscar el mejor (u óptimo) modelo : es decir, el que mejor se ajuste a los datos. Esto nos lleva a minimizar una función objetivo , es decir, un funcional que cuantifica cuán grandes son los residuos o cuán lejos están los datos predichos de los datos observados. Por supuesto, cuando tenemos datos perfectos (es decir, sin ruido), entonces el modelo recuperado debería ajustarse perfectamente a los datos observados. Una función objetivo estándar,${\displaystyle \varphi }$, tiene la forma: ${\displaystyle \varphi (p)=\|Fp-d_{\text{obs}}\|^{2}}$ dónde${\displaystyle \|\cdot \|}$es la norma euclidiana (será la${\displaystyle L^{2}}$norma cuando las mediciones son funciones en lugar de muestras) de los residuos. Este enfoque equivale a utilizar mínimos cuadrados ordinarios , un método ampliamente utilizado en estadística. Sin embargo, se sabe que la norma euclidiana es muy sensible a los valores atípicos: para evitar esta dificultad podemos pensar en utilizar otras distancias, por ejemplo la${\displaystyle L^{1}}$norma, en sustitución de la${\displaystyle L^{2}}$norma.

El enfoque probabilístico es muy similar al de mínimos cuadrados: si conocemos las estadísticas del ruido que contamina los datos, podemos buscar el modelo m más probable, que es el que cumple el criterio de máxima verosimilitud . Si el ruido es gaussiano , el criterio de máxima verosimilitud se convierte en un criterio de mínimos cuadrados, donde el producto escalar euclidiano en el espacio de datos se reemplaza por un producto escalar que involucra la covarianza del ruido. Además, si se dispone de información previa sobre los parámetros del modelo, se podría utilizar la inferencia bayesiana para formular la solución del problema inverso. Este enfoque se describe en detalle en el libro de Tarantola. ^{[ 10 ]} . El enfoque bayesiano para problemas inversos suele utilizar un proceso gaussiano para modelar la solución de las ecuaciones diferenciales (sin necesariamente linealizarlas) y emula (o entrena) el sistema mecanicista para que posteriormente actúe como una distribución a priori para la estimación (o calibración). ^{[ 11 ]}

Aquí utilizamos la norma euclidiana para cuantificar las discrepancias de los datos. Como se trata de un problema inverso lineal, la función objetivo es cuadrática. Para su minimización, es habitual calcular su gradiente utilizando el mismo razonamiento (como lo haríamos para minimizar una función de una sola variable). En el modelo óptimo${\displaystyle p_{\text{opt}}}$, este gradiente desaparece, lo que se puede escribir como: $${\displaystyle \nabla _{p}\varphi =2(F^{\mathrm {T} }Fp_{\text{opt}}-F^{\mathrm {T} }d_{\text{obs}})=0}$$ donde F ^T denota la transpuesta de la matriz F. Esta ecuación se simplifica a: $${\displaystyle F^{\mathrm {T} }Fp_{\text{opt}}=F^{\mathrm {T} }d_{\text{obs}}}$$

Esta expresión se conoce como la ecuación normal y nos da una posible solución al problema inverso. En nuestra matriz de ejemplo${\displaystyle F^{\mathrm {T} }F}$Resulta que generalmente tiene rango completo, por lo que la ecuación anterior tiene sentido y determina de forma única los parámetros del modelo: no necesitamos integrar información adicional para obtener una solución única.

Algunos problemas inversos tienen una solución muy simple, por ejemplo, cuando se tiene un conjunto de funciones unisolventes , es decir , un conjunto de${\displaystyle n}$⁠ funciones tales que evaluarlas en ⁠${\displaystyle n}$Los puntos distintos generan un conjunto devectores linealmente independientes . Esto significa que, dada una combinación lineal de estas funciones, los coeficientes se pueden calcular ordenando los vectores como columnas de una matriz y luego invirtiendo esta matriz. El ejemplo más simple de funciones unisolventes son los polinomios construidos, utilizando el teorema de unisolvencia , de manera que sean unisolventes. Concretamente, esto se hace invirtiendo la matriz de Vandermonde . Pero esta es una situación muy específica.

En general, la solución de un problema inverso requiere algoritmos de optimización sofisticados. Cuando el modelo se describe mediante un gran número de parámetros (el número de incógnitas involucradas en algunas aplicaciones de tomografía de difracción puede alcanzar los mil millones), resolver el sistema lineal asociado con las ecuaciones normales puede ser engorroso. El método numérico que se utilice para resolver el problema de optimización depende, en particular, del costo requerido para calcular la solución. ${\displaystyle Fp}$del problema directo. Una vez elegido el algoritmo apropiado para resolver el problema directo (una multiplicación matriz-vector directa puede no ser adecuada cuando la matriz${\displaystyle F}$es enorme), el algoritmo apropiado para llevar a cabo la minimización se puede encontrar en libros de texto que tratan sobre métodos numéricos para la solución de sistemas lineales y para la minimización de funciones cuadráticas (ver por ejemplo Ciarlet ^{[ 13 ]} o Nocedal ^{[ 14 ]} ).

Además, el usuario puede desear agregar restricciones físicas a los modelos: en este caso, debe estar familiarizado con los métodos de optimización con restricciones , un tema en sí mismo. En todos los casos, el cálculo del gradiente de la función objetivo suele ser un elemento clave para la solución del problema de optimización. Como se mencionó anteriormente, la información sobre la distribución espacial de un parámetro distribuido puede introducirse mediante la parametrización. También se puede considerar adaptar esta parametrización durante la optimización. ^{[ 15 ]}

Si la función objetivo se basa en una norma distinta de la norma euclidiana, debemos abandonar el ámbito de la optimización cuadrática. Como resultado, el problema de optimización se vuelve más difícil. En particular, cuando la${\displaystyle L^{1}}$La norma se utiliza para cuantificar el desajuste de los datos; la función objetivo ya no es diferenciable: su gradiente ya no tiene sentido. Entran en juego métodos específicos (véase, por ejemplo, Lemaréchal ^{[ 16 ] ) de optimización no diferenciable.}

Una vez calculado el modelo óptimo, debemos abordar la pregunta: "¿Podemos confiar en este modelo?" La pregunta se puede formular de la siguiente manera: ¿Qué tan grande es el conjunto de modelos que se ajustan a los datos "casi tan bien" como este modelo? En el caso de funciones objetivo cuadráticas, este conjunto está contenido en un hiperelipsoide, un subconjunto de${\displaystyle R^{M}}$(${\displaystyle M}$es el número de incógnitas), cuyo tamaño depende de lo que entendamos por "casi igual", es decir, del nivel de ruido. La dirección del eje más grande de este elipsoide ( vector propio asociado con el valor propio más pequeño de la matriz${\displaystyle F^{T}F}$) es la dirección de los componentes mal determinados: si seguimos esta dirección, podemos introducir una fuerte perturbación en el modelo sin cambiar significativamente el valor de la función objetivo y, por lo tanto, terminar con un modelo cuasi-óptimo significativamente diferente. Vemos claramente que la respuesta a la pregunta "¿podemos confiar en este modelo?" está gobernada por el nivel de ruido y por los valores propios del hessiano de la función objetivo o, equivalentemente, en el caso en que no se haya integrado ninguna regularización, por los valores singulares de la matriz${\displaystyle F}$Por supuesto, el uso de la regularización (u otro tipo de información previa) reduce el tamaño del conjunto de soluciones casi óptimas y, a su vez, aumenta la confianza que podemos depositar en la solución calculada.

Aquí nos centramos en la recuperación de un parámetro distribuido. Al buscar parámetros distribuidos, debemos discretizar estas funciones desconocidas. Al hacerlo, reducimos la dimensión del problema a un número finito. Pero ahora surge la pregunta: ¿existe alguna relación entre la solución que calculamos y la del problema inicial? Y otra pregunta: ¿qué entendemos por la solución del problema inicial? Dado que un número finito de datos no permite determinar un número infinito de incógnitas, la función de discrepancia de datos original debe regularizarse para garantizar la unicidad de la solución. Muchas veces, reducir las incógnitas a un espacio de dimensión finita proporciona una regularización adecuada: la solución calculada se asemejará a una versión discreta de la solución que buscábamos. Por ejemplo, una discretización simple suele funcionar para resolver el problema de deconvolución : funcionará siempre que no permitamos que aparezcan frecuencias faltantes en la solución numérica. Sin embargo, en muchos casos, la regularización debe integrarse explícitamente en la función objetivo.

Para comprender lo que puede suceder, debemos tener en cuenta que resolver un problema inverso lineal de este tipo equivale a resolver una ecuación integral de Fredholm de primera especie: $${\displaystyle d(x)=\int _{\Omega }K(x,y)p(y)dy}$$

dónde${\displaystyle K}$es el núcleo,${\displaystyle x}$y${\displaystyle y}$son vectores de${\displaystyle R^{2}}$, y${\displaystyle \Omega }$es un dominio en${\displaystyle R^{2}}$Esto se aplica a una aplicación 2D. Para una aplicación 3D, consideramos${\displaystyle x,y\in R^{3}}$. Tenga en cuenta que aquí los parámetros del modelo${\displaystyle p}$consiste en una función y que la respuesta de un modelo también consiste en una función denotada por${\displaystyle d(x)}$Esta ecuación es una extensión a dimensión infinita de la ecuación matricial.${\displaystyle d=Fp}$dado en el caso de problemas discretos.

Para que sea suficientemente suave${\displaystyle K}$El operador definido anteriormente es compacto en espacios de Banach razonables como el${\displaystyle L^{2}}$La teoría de F. Riesz establece que el conjunto de valores singulares de dicho operador contiene cero (de ahí la existencia de un espacio nulo), es finito o como máximo numerable y, en este último caso, constituyen una secuencia que tiende a cero. En el caso de un núcleo simétrico, tenemos una infinidad de autovalores y los autovectores asociados constituyen una base hilbertiana de${\displaystyle L^{2}}$Por lo tanto, cualquier solución de esta ecuación está determinada salvo una función aditiva en el espacio nulo y, en el caso de infinitos valores singulares, la solución (que involucra el recíproco de autovalores arbitrariamente pequeños) es inestable: ¡dos ingredientes que hacen de la solución de esta ecuación integral un típico problema mal condicionado! Sin embargo, podemos definir una solución a través de la pseudoinversa del mapeo directo (nuevamente salvo una función aditiva arbitraria). Cuando el mapeo directo es compacto, la regularización clásica de Tikhonov funcionará si la usamos para integrar información previa que establece que${\displaystyle L^{2}}$La norma de la solución debe ser lo más pequeña posible: esto hará que el problema inverso esté bien planteado. Sin embargo, como en el caso de dimensión finita, debemos cuestionar la confianza que podemos depositar en la solución calculada. Nuevamente, básicamente, la información reside en los autovalores del operador hessiano. Si se exploran subespacios que contienen autovectores asociados a autovalores pequeños para calcular la solución, entonces difícilmente se podrá confiar en ella: algunos de sus componentes estarán mal determinados. El autovalor más pequeño es igual al peso introducido en la regularización de Tikhonov.

Los núcleos irregulares pueden producir un mapa hacia adelante que no es compacto e incluso no está acotado si ingenuamente equipamos el espacio de modelos con el${\displaystyle L^{2}}$norma. En tales casos, el hessiano no es un operador acotado y la noción de autovalor deja de tener sentido. Se requiere un análisis matemático para convertirlo en un operador acotado y diseñar un problema bien planteado: una ilustración se puede encontrar en ^{[ 17 ]} . Nuevamente, debemos cuestionar la confianza que podemos depositar en la solución calculada y debemos generalizar la noción de autovalor para obtener la respuesta. ^{[ 18 ]}

El análisis del espectro del operador hessiano es, por lo tanto, un elemento clave para determinar la fiabilidad de la solución calculada. Sin embargo, dicho análisis suele ser una tarea muy compleja. Esto ha llevado a varios autores a investigar enfoques alternativos en el caso de que no nos interesen todos los componentes de la función desconocida, sino solo los sub-desconocidos que son las imágenes de la función desconocida mediante un operador lineal. Estos enfoques se conocen como el método de Backus y Gilbert ^{[ 19 ]} , el método de los centinelas de Lions ^{[ 20 ]} y el método SOLA ^{[ 21 ]} . Estos enfoques resultaron estar fuertemente relacionados entre sí, como se explica en Chavent ^{[ 22 ].} Finalmente, el concepto de resolución limitada , frecuentemente invocado por los físicos, no es más que una visión específica del hecho de que algunos componentes mal determinados pueden corromper la solución. Pero, en términos generales, estos componentes mal determinados del modelo no están necesariamente asociados con altas frecuencias.

El objetivo de la deconvolución es reconstruir la imagen o señal original.${\displaystyle p(x)}$que aparece como ruidoso y borroso en los datos.${\displaystyle d(x)}$. ^{[ 23 ]} Desde un punto de vista matemático, el núcleo${\displaystyle K(x,y)}$Aquí solo depende de la diferencia entre${\displaystyle x}$y${\displaystyle y}$.

En estos métodos intentamos recuperar un parámetro distribuido, cuya observación consiste en la medición de las integrales de este parámetro realizadas a lo largo de una familia de líneas. Denotamos por${\displaystyle \Gamma _{x}}$la línea en esta familia asociada con el punto de medición${\displaystyle x}$La observación en${\displaystyle x}$Por lo tanto, se puede escribir como: $${\displaystyle d(x)=\int _{\Gamma _{x}}w(x,y)p(y)\,dy}$$ dónde${\displaystyle w(x,y)}$una función de ponderación conocida. Comparando esta ecuación con la integral de Fredholm anterior, observamos que el núcleo${\displaystyle K(x,y)}$es una especie de función delta que alcanza su máximo en la línea${\displaystyle {\Gamma _{x}}}$Con un núcleo de este tipo, el mapa hacia adelante no es compacto.

Las primeras publicaciones sobre conducción de calor inversa surgieron de la determinación del flujo de calor superficial durante la reentrada atmosférica a partir de sensores de temperatura enterrados. ^{[ 33 ] [ 34 ]} Otras aplicaciones donde se necesita el flujo de calor superficial pero los sensores superficiales no son prácticos incluyen: dentro de motores alternativos, dentro de motores de cohetes; y, pruebas de componentes de reactores nucleares. ^{[ 35 ]} Se han desarrollado diversas técnicas numéricas para abordar la mala condición y la sensibilidad al error de medición causado por la amortiguación y el retardo en la señal de temperatura. ^{[ 36 ] [ 37 ] [ 38 ]}

Si bien los problemas inversos lineales se resolvieron completamente desde el punto de vista teórico a finales del siglo XIX , solo una clase de problemas inversos no lineales lo estaba antes de 1970: los problemas inversos espectrales y de dispersión (en una dimensión espacial) , tras la obra fundamental de la escuela matemática rusa ( Krein , Gelfand , Levitan, Marchenko ). Chadan y Sabatier realizaron una amplia revisión de los resultados en su libro "Problemas inversos de la teoría de la dispersión cuántica" (dos ediciones en inglés, una en ruso).

En este tipo de problemas, los datos son propiedades del espectro de un operador lineal que describe la dispersión. El espectro está compuesto por autovalores y autofunciones , que forman el "espectro discreto", y generalizaciones, denominadas espectro continuo. Lo más destacable es que los experimentos de dispersión solo proporcionan información sobre el espectro continuo, y que conocer su espectro completo es necesario y suficiente para recuperar el operador de dispersión. Por lo tanto, tenemos parámetros invisibles, mucho más interesantes que el espacio nulo, que posee una propiedad similar en problemas inversos lineales. Además, existen movimientos físicos en los que el espectro de dicho operador se conserva como consecuencia del movimiento. Este fenómeno se rige por ecuaciones diferenciales parciales no lineales especiales, como la ecuación de Korteweg-de Vries . Si el espectro del operador se reduce a un único autovalor, su movimiento correspondiente es el de una onda solitaria que se propaga a velocidad constante y sin deformación, denominada " solitón ".

La obtención de una señal perfecta y sus generalizaciones para la ecuación de Korteweg-de Vries u otras ecuaciones diferenciales parciales no lineales integrables son de gran interés, con numerosas aplicaciones potenciales. Esta área se estudia como una rama de la física matemática desde la década de 1970. Los problemas inversos no lineales también se estudian actualmente en muchos campos de la ciencia aplicada (acústica, mecánica, mecánica cuántica, dispersión electromagnética, en particular sondeos de radar, sondeos sísmicos y prácticamente todas las modalidades de imagen).

Un último ejemplo relacionado con la hipótesis de Riemann fue dado por Wu y Sprung, la idea es que en la antigua teoría cuántica semiclasica el inverso del potencial dentro del hamiltoniano es proporcional a la semiderivada de la función de conteo de autovalores (energías) n ( x ). 

El objetivo es recuperar el coeficiente de difusión en la ecuación diferencial parcial parabólica que modela flujos de fluidos monofásicos en medios porosos. Este problema ha sido objeto de numerosos estudios desde un trabajo pionero realizado a principios de los años setenta. ^{[ 39 ]} En cuanto a los flujos bifásicos, un problema importante es estimar las permeabilidades relativas y las presiones capilares. ^{[ 40 ]}

El objetivo es recuperar las velocidades de las ondas (ondas P y S) y las distribuciones de densidad a partir de los sismogramas . Estos problemas inversos son de gran interés en sismología y geofísica de exploración . Básicamente, podemos considerar dos modelos matemáticos:

• La ecuación de onda acústica (en la que se ignoran las ondas S cuando las dimensiones espaciales son 2 o 3)
• La ecuación de elastodinámica en la que las velocidades de las ondas P y S se pueden derivar de los parámetros de Lamé y la densidad.

Estas ecuaciones hiperbólicas básicas pueden mejorarse incorporando atenuación , anisotropía y otros efectos físicos.

La solución del problema inverso en la ecuación de onda 1D ha sido objeto de muchos estudios. Es uno de los pocos problemas inversos no lineales para los que podemos demostrar la unicidad de la solución. ^{[ 8 ]} El análisis de la estabilidad de la solución fue otro desafío. ^{[ 41 ]} Se desarrollaron aplicaciones prácticas, utilizando el enfoque de mínimos cuadrados. ^{[ 41 ] [ 42 ]} Se intentó la extensión a problemas 2D o 3D y a las ecuaciones de elastodinámica desde los años 80, ¡pero resultó ser muy difícil  ! Este problema, a menudo denominado Inversión de Forma de Onda Completa (FWI), aún no está completamente resuelto: entre las principales dificultades se encuentran la existencia de ruido no gaussiano en los sismogramas, problemas de salto de ciclo (también conocido como ambigüedad de fase) y el comportamiento caótico de la función de error de datos. ^{[ 43 ]} Algunos autores han investigado la posibilidad de reformular el problema inverso para que la función objetivo sea menos caótica que la función de error de datos. ^{[ 44 ] [ 45 ]}

Al comprender la dificultad del problema inverso en la ecuación de onda, los sismólogos investigaron un enfoque simplificado que utiliza la óptica geométrica. En concreto, se propusieron invertir la distribución de velocidad de propagación, conociendo los tiempos de llegada de los frentes de onda observados en los sismogramas. Estos frentes de onda pueden asociarse con llegadas directas o con reflexiones asociadas a reflectores cuya geometría debe determinarse junto con la distribución de velocidad.

La distribución del tiempo de llegada${\displaystyle {\tau }(x)}$(${\displaystyle x}$es un punto en el espacio físico) de un frente de onda emitido desde una fuente puntual, satisface la ecuación de Eikonal : $${\displaystyle \|\nabla \tau (x)\|=s(x),}$$ dónde${\displaystyle s(x)}$denota la distribución de lentitud (recíproca de la velocidad). La presencia de${\displaystyle \|\cdot \|}$Esto hace que esta ecuación sea no lineal. Clásicamente se resuelve mediante el lanzamiento de rayos (trayectorias respecto de las cuales el tiempo de llegada es constante) desde la fuente puntual.

Este problema es similar a la tomografía: los tiempos de llegada medidos son la integral de la lentitud a lo largo de la trayectoria del rayo. Sin embargo, este problema, similar a la tomografía, es no lineal, principalmente porque la geometría desconocida de la trayectoria del rayo depende de la distribución de velocidad (o lentitud). A pesar de su carácter no lineal, la tomografía de tiempos de viaje ha demostrado ser muy eficaz para determinar la velocidad de propagación en la Tierra o en el subsuelo, siendo este último aspecto un elemento clave para la obtención de imágenes sísmicas, en particular utilizando los métodos mencionados en la sección "Tomografía de difracción".

Dado que el mapeo directo es no lineal, es probable que la función de error de datos no sea convexa, lo que hace que las técnicas de minimización local sean ineficientes. Se han investigado varios enfoques para superar esta dificultad:

• uso de técnicas de optimización global como el muestreo de la función de densidad posterior y el algoritmo de Metropolis en el marco probabilístico del problema inverso, ^{[ 47 ]} algoritmos genéticos (solos o en combinación con el algoritmo de Metropolis: ver ^{[ 48 ]} para una aplicación a la determinación de permeabilidades que coinciden con los datos de permeabilidad existentes), redes neuronales, técnicas de regularización incluido el análisis multiescala;
• reformulación de la función objetivo de mínimos cuadrados para hacerla más suave (ver ^{[ 44 ] [ 45 ]} para el problema inverso en las ecuaciones de onda).

Los problemas inversos, especialmente en dimensión infinita, pueden ser de gran tamaño, lo que requiere un tiempo de cálculo considerable. Cuando el mapeo directo es no lineal, las dificultades computacionales aumentan y minimizar la función objetivo puede resultar complicado. A diferencia de la situación lineal, el uso explícito de la matriz hessiana para resolver las ecuaciones normales no tiene sentido en este caso, ya que la matriz hessiana varía según el modelo. Para algunos modelos, resulta mucho más eficaz la evaluación del gradiente de la función objetivo. Se puede ahorrar un esfuerzo computacional importante evitando el cálculo muy pesado del jacobiano (a menudo denominado " derivadas de Fréchet "): el método del estado adjunto, propuesto por Chavent y Lions, ^{[ 49 ]} busca evitar este cálculo tan pesado. Actualmente se utiliza ampliamente. ^{[ 50 ]}