TUTORIAL UNA REVISIÓN DE INVERSIÓN POR MÍNIMOS CUADRADOS Y SU APLICACIÓN

PARA PROBLEMAS GEOFÍSICOS

L.R. LINES y S. TREITEL
Abstract
La inversión geofísica involucra la estimación de los parámetros de un modelo de la tierra
postulado de un conjunto de observaciones. Dado que la respuesta del modelo asociado
puede ser funciones no lineales de los parámetros del modelo, la técnica de mínimos
cuadrados no lineales demuestra ser útil para representar la inversión. Un caso común de
inversión aplica mínimos cuadrados de amortiguadores iterativos a través del uso del
método de Marquardt-Levenberg. Tradicionalmente, este método ha sido implementado
para resolver la ecuación normal asociada en procedimientos convencionales. Sin embargo,
la Descomposición del Valor Singular (SVD por sus siglas en inglés) produces significantes
mejoras en precisión computacional cuando aplica para el mismo sistema de ecuaciones
normal. La modelación de mínimos cuadrados iterativos tiene aplicación en una amplia
variedad de problemas geofísicos. Los siguientes dos ejemplos ilustran el enfoque: (1)
deconvolución de la onda sísmica, y (2) la localización de una cuña enterrada a partir de
datos de gravimetría en superficie. Mas generalmente, la inversión de mínimos cuadrados
puede ser usado para estimar modelos de la tierra para un conjunto de observaciones
geofísicas para los cual una descripción matemática apropiada está disponible.
1. Problemas de inversión geofísica
La inversión geofísica podría ser vista como un intento para ajustar la respuesta de un
modelo de la tierra de superficie idealizada para un conjunto finito de observaciones
actuales. Distinguimos entre un modelo, sus parámetros del modelo, y la respuesta del
modelo asociado. Un modelo consiste de un conjunto de relaciones representando una
abstracción matemática particular de un proceso observado. Esas ecuaciones en turno
dependen de un cierto número, llamado p, de parámetros modelo los cuales requerimos
para estimar a partir de los datos actuales. La respuesta del modelo consiste de los datos
sintéticos producidos por una realización particular del modelo. El propósito del modelo es
extraer los parámetros estimados de un intento de ajuste del modelo respuesta a los datos
observados. Por ejemplo, el medio acústico estratificado de una dimensión obedece a la
ecuación de la onda acústica siendo un modelo, cuyos parámetros del modelo son las
impedancias y espesores de las capas, y cuyo modelo de respuesta es el familiar sismograma
sintético impulsivo de incidencia normal. En el riesgo de exceder lo obvio, señalamos que si
el modelo es físicamente irrelevante, los parámetros estimados que podríamos extraer de
ello son insignificantes. Puesto que los datos geofísicos actuales son de alcance finito y
porque ellos siempre contienen componentes los cuales no podemos explicar en términos
geológicos, la inversión puede nunca ser única, eso es, más de una solución satisfará las
observaciones dentro de un erro prescrito.
En este tratamiento no hacemos pretender a un estudio exhaustivo del arte de la inversión
geofísica. Más bien, hemos elegido el método de mínimos cuadrados debido a su robustez
matemática cuando los datos son registrados, en las palabras de Jackson (1972), “inexacto,
insuficiente e inconsistente”. Numerosos estudios pertinentes han aparecido durante las
últimas dos décadas. Ello incluye la inversión de datos de la tierra global para uso libre de
oscilaciones (Crosson 1976, Neumann, Gjoystdal and Ursin 1981), la inversión de datos
electromagnéticos (Inman 1975, Jupp and Vozoff 1975, Oristaglio and Worthington 1980),
and la inversión de datos de gravedad (Oldenburg 1974, Vigneresse 1977). Un compendio
integral de métodos inversos para sensores remotos ha sido dado por Twomey (1977),
mientras que excelentes tratamientos de los conceptos matemáticos subyacentes han sido
publicados por Jackson (1972, 1979) y Aki y Richards (1980).
Quizás la respuesta de los modelos geofísicos son funciones no lineales de los parámetros
del modelo apropiado. El también llamado método de “Marquardt-Levenberg” nos provee
con un poderoso y versátil método para estimas tales parámetros en una manera iterativa
(Levenberg, 1944, Marquardt, 1963). Nuestro artículo es un intento para describir los
aspectos notables de la inversión por mínimos cuadrados en términos relativamente
simples. Sin embargo, no podemos evitar el uso de algebra vectorial y matricial, asi que
debemos asumir que el lector tiene la capacidad para sobrellevarnos.
2. Inversión por mínimos cuadrados
Una respuesta del modelo puede ser tanto una función lineal como una no lineal de los
parámetros del modelo. Entonces para el modelo de convolución lineal yt=at*bt, donde *
denota convolución, el dato yt es una función lineal de los parámetros del modelo at y bt,
pero si, suponemos que, yt=at2, entonces yt es una función lineal de at. Una discusión
compresiva de la relación entre métodos de inversión por mínimos cuadrados lineal y no
lineal pueden ser encontrados en el libro de Draper y Smith (1981, capitulo 10). Un
comentario sobre la terminología está en orden: Draper y Smith son estadísticos y usan la
palabra “regresión” cuando nosotros usamos “inversión”. Para complicar más allá el asunto,
los ingenieros eléctricos frecuentemente refieren a esos métodos como “identificación de
sistemas” o “identificación del parámetro”. Otros términos están también en uso.
2.1. Inversión por mínimos cuadrados lineal
Antes de proceder a los problemas de inversión no lineal, examinaremos el problema
inverso lineal sencillo. La extensión de esas ideas al caso no lineal es introducida en la
sección 2.3.
La estrategia básica es minimizar la suma de cuadrados de los errores entre la respuesta del
modelo y las observaciones. Las n observaciones para un conjunto de datos de la tierra es
representado por el vector
y= col (y1, y2,…, yn),
y la respuesta del modelo es el vector
f= col (f1, f2,…, fn).
El modelo es una función de p parámetros los cuales son elementos de un vector Ѳ,
Ѳ = col (Ѳ 1, Ѳ 2,…, Ѳ p).
Ѳj0 es un estimado de el parámetro Ѳj, (j=1,…, p), y f0 es la respuesta del modelo inicial. Si
la respuesta del modelo f es una función lineal de los parámetros, una perturbación de la
respuesta del modelo respecto a Ѳ0 puede ser representado por la expansión de Taylor de
primer orden:

o, en notación matricial,
f= f0 + Zẟ, (5)
donde Z es la matriz Jacobiana n x p de derivadas parciales con elementos

y ẟ= Ѳ - Ѳj0 es el parámetro del vector de cambio con elementos ẟj representando los

cambios, o perturbaciones en los parámetros Ѳj, eso es
ẟj= Ѳj - Ѳj0, (j=1,…, p).
Nuestra opción de perturbaciones en Ѳ estará hecha para que minimice la suma de los
cuadrados de los errores entre la respuesta del modelo y los datos. e representa el vector
de error expresando la diferencia entre el modelo respuesta f y los datos observados y:
y – f= e. (7)
Combinando (5) y (7) se produce
y – (f0 + Zẟ)= e
o
y – f0= Zẟ + e. (8)
El vector y – f0, el cual contiene las diferencias entre la respuesta modelo inicial y los datos
observados, es llamado el vector de discrepancia g, tal que
g= y – f0
y
e= g - Zẟ. (9)
Antes describimos un método de mínimos cuadrados dirigido a un estimado de los
parámetros del vector de cambio ẟ notamos que los problemas inversos geofísicos están
generalmente no bien planteados, esos es, que la matriz Jacobiana Z de n x p no es
usualmente cuadrado y de rango completo. De hecho, mucho de esos problemas son
sobredeterminados, eso es, el número de puntos de dato excede el número de parámetros
del modelo, i.e., n > p; la matriz Z podría o no ser de rango completo. En el caso de rango
deficiente, existen técnicas para mejorar el problema resultante. El método de Marquardt-
Levenberg (ver sección 2.2) puede lograr esto muy bien. Un trato claro del problema ha sido
dado por Lanczos (1961, capitulo III).
En el más simple enfoque de mínimos cuadrados o “Gauss-Newton”, pretendimos
minimizar el error cuadrado acumulativo S= eT e con respecto al parámetro de vector de
cambio ẟ.
De (9), tenemos
S= eT e= (g - Zẟ)T (g - Zẟ). (10)
Minimización de S con respecto a ẟ requiere que
𝜕𝑆
= 0. (11)
𝜕ẟ

(La diferenciación con respecto a un vector es descrito por Graybill (1969) e implica que
𝜕𝑆/𝜕ẟi = 0 para toda i.)
Sustituyendo (10) en (11) da
𝜕
(ẟTZTZẟ - gTZẟ - ẟTZTg + gTg)= 0.
𝜕ẟ

Llevando a cabo la diferenciación con respecto a ẟ, obtenemos las también llamadas

“ecuaciones normales”
ZTZẟ = ZTg,
cuya solución para los parámetros del vector de cambio ẟ es
ẟ= (ZTZ)-1 ZTg, (14)
la cual es también conocida como la solución de Gauss-Newton. Este conjunto de
ecuaciones ha tenido muchas importantes aplicaciones en la teoría de inversión y filtros
digitales. Los resultados anteriores pueden ser derivados a partir de netamente
consideraciones geométricas, ver, por ejemplo, Strang (1980, capitulo 3). De hecho, el
término “normal” surge de la propiedad geométrica que el vector e de error de los mínimos
cuadrados es perpendicular, o normal a la columna de vectores de la matriz Z.

El valor actual del error de mínimos cuadrados acumulativos 𝑆̂ resulta de sustituir la

solución de mínimos cuadrados (14) en (10),

𝑆̂= (Zẟ̂ - g)T(Zẟ̂ – g)

donde ẟ̂ es la solución de los mínimos cuadrados (14). Sustituyendo esta expresión para ẟ̂
en la ecuación anterior y llevando a cabo algún algebra matricial se obtiene

𝑆̂= gT(In – ZZL-1) g,

donde ZL-1= (ZTZ)-1 ZT. La matriz ZL-1 de p x n es el mínimo cuadrado, o la inversión Lanczos
de Z (ver también sección 3.1), y In es la matriz identidad de n x n. Observamos que 𝑆̂ llega
a ser pequeño en cualquier momento que el producto ZZL-1 aproxima a In. En la literatura la
inversión Lanczos es también frecuentemente referida como la “natural” o “inversión
generalizada” de la matriz rectangular Z de n x p.
La solución de mínimos cuadrados (14) ha resultado de la directa minimización del error
cuadrado acumulativo S. Este también llamado solución de mínimos cuadrados
“espontáneo” debido a que tiene algunas propiedades indeseables, las cuales
describiremos.
Una dificultad obvia ocurre cuando la inversa de ZTZ no existe, eso es, cuando la matriz ZTZ
es singular. Aun si ZTZ existe, podemos enfrentarnos con una solución divergente, o lidiar
con convergencia baja. Esto puede suceder en cualquier momento que la estimación inicial
del modelo, f0, sea pobre (Draper y Smith 1981). Tan pronto como (ZTZ) llega a ser
cercanamente singular, los elementos del vector solución ẟ tiende a crecer sin límite. Smith
y Shanno (1971) han descrito como una ZTZ casi singular tiende a catapultar el vector de
parámetros actualizados Ѳ= Ѳ0 + ẟ lejos de una solución aceptable.
2.2. El método de Marquardt-Levenberg
Para reducir las dificultades cuando la matriz ZTZ es casi singular, resolvemos un problema
de mínimos cuadrados alternado. Imponemos la condición obligada de que la suma de los
cuadrados, o energía de los elementos de los parámetros del vector de cambio ẟ, sea
definido para una cantidad finita, es decir ẟ02. Este enfoque fue introducido por Levenberg
(1944) y más tarde descrito a detalle por Marquardt (1963). Es conocido algunas veces como
el método de mínimos amortiguados; otros refieren a esto como “regresión de cresta”
(Inman 1975). El efecto de este contraste es impedir oscilaciones ilimitadas en la solución;
i.e., suavizar los parámetros del vector de cambio ẟ.
La solución de mínimos cuadrados forzada surge para resolver un problema de los
multiplicadores de Lagrange en el cual eTe es el argumento minimizado al forzar que ẟT ẟ=
ẟ02. Entonces, elegimos ẟ para minimizar una función de coste S (ẟ,β),
S (ẟ, β)= eTe + β (ẟT ẟ - ẟ02),
donde β es un multiplicador de Lagrange.
La diferenciación con respecto al vector ẟ surge de una forma modificada de ecuación
normal
(ZTZ + βI) ẟ= ZTg,
asi que
ẟ= (ZTZ + βI)-1 ZTg, (16b)
La comparación de (14) con (16b) muestra que la restricción ha producido un método para
evitar singularidades o casi singularidades en la matriz ZTZ. Añadiendo una constante, β, a
la diagonal principal de ZTZ, tenemos efectivamente agregada un nivel d.c. a los
eigenvalores de la matriz ZTZ así que ninguno de los eigenvalores puede desaparecer.
Levenberg (1944) llama al multiplicador de Lagrange β un “factor de amortiguamiento”, ya
que efectivamente amortigua los cambios en el vector de parámetros Ѳ para limitar la
energía en los parámetros del vector de discrepancia ẟ.
La solución (16b) tiene otras varias características interesantes. Esto es hibrido, porque
combina el también llamado “método de steepest descent” (método del descenso más
empinado) con el método de mínimos cuadrados. La solución de steepest descent es normal
para un perfil de función de coste dado, para el cual S (Ѳ)= eTe= constante, o dS (Ѳ)= 0
(Smith y Shanno 1971). Esta condición es satisfecha por el vector de columna steepest
descent con p componentes
ẟg= -𝛁 S (Ѳ),
donde 𝛁 es el operador de gradiente 𝜕/𝜕Ѳj, j =1,…, p. Notamos que +𝛁 S (Ѳ) es entonces
el vector steepest ascent.
Entonces

el elemento jth del vector columna 𝛁 S (Ѳ) es

Debido a que ei= yi –fi (ver (7)), tenemos

y por lo tanto el vector columna ẟg= -𝛁 S (Ѳ) es
ẟg= 2ZTe
o
ẟg= 2ZT (y - f),
donde ZT es la transpuesta de la matriz Z Jacobiana de n x p. La iteración comienza con un
modelo inicial cuya respuesta es f= f0, así
ẟg= 2ZT (y – f0),
pero y – f0= g (ver (8)), y el vector steepest descent es
ẟg= 2ZTg. (18)
Dado que ẟg esta en la dirección de decrecer S, la convergencia tiende a ocurrir, pero
usualmente en un bajo rango. Además, problemas computacionales pueden ocurrir cuando
el tamaño del grado ẟj= Ѳj – Ѳj0 llega a ser también pequeño.
Generalmente hablando, el método steepest descent es óptimo cuando S (Ѳ) es grande,
mientras que el método de mínimos cuadrados llega a ser efectivo cuando S (Ѳ) es pequeño.
El enfoque de Marquardt-Levenberg entonces aprovecha el uso de ambos métodos. La
técnica determina los parámetros del vector de cambio ẟ en cada grado en la iteración.
Marquardt (1963) demostró que los vectores ẟ y ẟg deben de estar dentro de un máximo
de 90° uno de otro, y esa es una propiedad necesaria para converger. La comparación de
(16b) con (18) muestra como β llega a ser grande ẟ= β-1 ZTg, i.e., ẟ llega a ser proporcional
al vector gradiente ẟ. Sin embargo, como β aumenta y aumenta, los elementos de los
parámetros del vector de cambio ẟ disminuyen y disminuyen. Marquardt también mostro
que el ángulo entre ẟ y ẟg es monotónicamente un función decreciente de β.
Esos aspectos son mostrados en la figura 1 para un modelo con dos parámetros Ѳ1 y Ѳ2. En
un caso tan simple las curvas de constante error cuadrado acumulativo S (Ѳ 1 ,Ѳ2) surgen
como contornos en el plano Ѳ1 - Ѳ2. Suponiendo que nuestro punto de inicio es
Ѳ0= (Ѳ10 ,Ѳ20),
y suponiendo que el problema tiene solo un mínimo único en el punto
̂ = (Ѳ1
Ѳ ̂,̂Ѳ2),
el es entonces la solución de los mínimos cuadrados iterativos deseada. La figura 1 ilustra
el primer grado en nuestro intento por alcanzar este mínimo cuando usamos tanto el
steepest descent, el Marquardt-Levenberg, o la aproximación de Gauss-Newton.
Observemos que el vector steepest descent es perpendicular al error más alejado del
contorno en el punto Ѳ0= (Ѳ10 ,Ѳ20). El vector Marquardt-Levenberg es asociado por los
vectores Gauss-Newton y steepest descent, y de una manera ofrece un compromiso entre
cada extremo. La aplicación sucesiva de alguno de esos algoritmos eventualmente alcanzara

el mínimo en (Ѳ1̂ ,̂Ѳ2); en general, el método de Marquardt-Levenberg lo hará más

eficientemente (ver también sección 4.2).