Ejemplo sencillo para ajustar polinomios de tendencia 2D

Métodos potenciales de prospección, FCAG, 2024.

Lectura del dato

Vamos a utilizar el archivo datos2D.dat de anomalías de Bouguer dispuestos en una cuadrícula. En primer lugar utilizamos algunos comandos GNU/Linux en la terminal para observar la estructura del dato.

$> head "datos2D.dat"

# x[km]      y[km]    AB[mGal]

 0	0	8.36
 1	0	8.46
 2	0	8.51
 3	0	8.56
 4	0	8.72
 5	0	8.67
 6	0	8.79
 7	0	8.78

Viendo el encabezado, notamos que el dato está presentado en la forma

$x$ (km)	$y$ (km)	$g_z$ (mGal)
...	...	...
...	...	...
...	...	...

Ahora procedemos a leer el archivo para procesar. Para ello podemos escribir:

data     = np.loadtxt('datos2D.dat')       # las líneas que comienzan con '#' en el archivo son ignoradas
x, y, gz = data[:,0], data[:,1], data[:,2] # columnas del dato en cada vector

N  = len(gz)           # número de observaciones de AAB [mgal]
nx = len(np.unique(x)) # número de coordenadas únicas en x
ny = len(np.unique(y)) # número de coordenadas únicas en y

print(f"Contamos con {N} observaciones.")
print(f"nx={nx}, ny={ny}, (nx * ny = {nx*ny}).")

GZ = gz.reshape(ny,nx) # genero matriz con el vector de mediciones de gravedad (sólo para graficar)

Contamos con 176 observaciones.
nx=16, ny=11, (nx * ny = 176).

Ajuste de polinomios de tendencia

Ajustamos al dato $g_z$ polinomios de grado $M$,

\begin{equation} \hat{g}_z(x,y;M) = \begin{cases} c, & \; M=0;\\ ax + by + c, & \; M=1;\\ ax^2 + by^2 + cxy + dx + ey + f, & \; M=2. \end{cases} \end{equation}

Calculamos la matriz de diseño $\mathbf{A}$ para $M=1$:

x    = x.reshape(N,1)
y    = y.reshape(N,1)
ones = np.ones((N,1))
A    = np.hstack((x, y, ones))

print(f"Tamaño de A: {np.shape(A)}")

Tamaño de A: (176, 3)

Resolvemos el sistema $\mathbf{A}\mathbf{m}=\mathbf{g^o}$, para $\mathbf{m}$ por cuadrados mínimos. Aquí $\mathbf{g^o}$ es el dato observado dispuesto en un vector columna. La solución dada por cuadrados mínimos la denotamos $\mathbf{m}^*$ y en código m_estimated:

from scipy.linalg import lstsq
m_estimated,_,_,_ = lstsq(A,gz)

Grado del polinomio de ajuste, M = 1
Vector de coeficientes, m*       = [-0.022066 -0.209915 8.734554]

Como valor agregado, podemos verificar el resultado si utilizamos un algoritmo de gradientes conjugados (CG) programado en una clase preliminar:

def CG(A,b):
    ...
    return m

m_estimated = CG(A.T @ A,A.T @ gz)

CG: Según la tolerancia 1e-10, convergió en la iteración 4
vector de coeficientes, m* =  [-0.022066 -0.209915 8.734554]

Calculo de la superficie de tendencia y el mapa residual

Con el vector de coeficientes $\mathbf{m^*}$, obtenemos por medio del problema directo la superficie de tendencia en los puntos observados:

$$\mathbf{\hat{g}} = \mathbf{A}\mathbf{m^*},$$

y el campo residual,

$$\Delta\mathbf{g} = \mathbf{g^o}-\mathbf{\hat{g}}.$$

En python podemos generar la superficie de tendencia y el campo residual así:

# Superficie de tendencia:
GZ_fit = (A @ m) # en python 3, '@' es producto matricial
GZ_fit = GZ_fit.reshape(ny,nx) # construyo cuadrícula del ajuste 
# Campo residual (resta de cuadrículas):
diff = GZ - GZ_fit

Resultados

Visualizamos el dato original $g_z$, la superficie ajustada según mínimos cuadrados $\hat{g}_z$ y el residuo $\Delta g_z=g_z-\hat{g}_z$. Observamos que las líneas de contorno del mapa residual no son suaves debido a que el mismo no ha sido interpolado. En los gráficos, utilizamos la misma escala vertical para el dato original y la superficie de tendencia. Para visualizar el campo residual, utilizamos una escala diferente centrada en 0 mGal.

Eso es todo por hoy.