Algoritmo para la estimación de la turbidez sobre el Río Paraná

Autor/a

Afiliación

Víctor Gustavo Gómez

GISTAQ (UTN-FRRe)

Resumen

Este sitio web contiene información sobre la estimación de la turbidez por teledetección en la cuenca media del Río Paraná. La turbidez es uno parámetros de interés dentro proyecto Estimar indicadores de calidad de agua en la cuenca media del río Paraná para el desarrollo de un algoritmo mediante técnicas de teledetección satelital (MSECRE0008604), desarrollado por el Grupo de Investigación Sobre Temas Ambientales y Químicos (GISTAQ) de la Universidad Tecnológica Nacional Facultad Regional Resistencia (UTN-FRRe).

Se utilizarán imágenes del satélite Sentinel-2 con corrección automática, de las cuales se obtiene la reflectancia de superficie del agua. Se buscará la relación entre la reflectancia y la turbidez por métodos de regresión tradionales y machine learning. Una vez obtenido el algoritmo que relacione ambas propiedades, se desarrollaran mapas de distribución espacial.

Palabras clave

GISTAQ, UTN, FRRe, Algoritmo, Turbidez, Machine learning, Teledetección

1 Turbidez

La turbidez se refiere a la opacidad o falta de claridad en un líquido provocada por la presencia de partículas suspendidas. Este fenómeno es un indicador clave en el monitoreo de la calidad del agua y su influencia en diferentes ecosistemas es significativa.

La turbidez es un indicador de la calidad del agua, reflejando la presencia de partículas en suspensión. Su medición es crucial para garantizar la potabilidad del agua y la salud de los ecosistemas acuáticos. Este fenómeno puede ser resultado de diversas causas, como la erosión del suelo, la actividad biológica y la contaminación. La comprensión de la turbidez y su impacto es esencial para la gestión de recursos hídricos y la protección del medio ambiente.

La turbidez viene determinada por la dispersión de la luz causada por la materia suspendida en el agua, se obtiene normalmente mediante un turbidímetro, que proporciona medidas en Nephelometric Turbidity Unit (NTU) y mide la dispersión de un rayo de luz en el agua a 90º de la luz incidente [1].

Muchas propiedades, como la clorofila-a (Chl-a), sólidos suspendidos totales (SST) y la materia orgánica disuelta coloreada (CDOM), se utilizan a menudo como indicadores del estado del agua. Estos constituyentes del agua a su vez son responsables de la turbidez.

Existe una fuerte correlación entre turbidez y sólidos suspendidos totales, por lo que se puede estimar SST a partir de la turbidez. Por lo general, es una relación directa, a mayor concentración de SST mayor turbidez.

Existe una relación inversa entre la Turbidez y la profundidad del disco de Secchi (a valores bajos de secchi mayor turbidez), por lo que también se puede estimar turbidez a partir de mediciones de disco de secchi.

1.1 Métodos tradicionales

Tabla 1: Características principales de algoritmos tradicionales para la estimación de turbidez.

Ecuación	Bandas (nm)	Métricas	Aguas	Plataforma	Referencia
\(1.559e^{35.533 \cdot B03} \\ 1.879e^{37.745(B03 \cdot B5)/(B04+B12)}\)	B03, B04, B05, B12	\(R^{2}\), RMSE, MAE	Lago1	Sentinel-2	[2]
\(2677.2 \cdot B04^{1.856}\)	B04	\(R^{2}\), RMSE, Bias	Interiores variadas2	Landsat-8	[3]
\(969-1.5468 \cdot R_{1200nm}+2.07 \frac{B8A}{B02}\)	B02, B8A, 1200nm	IOA, SI, RMSE, MAE	Río3	Landsat-8	[4]
\(y=-1.1+5.8 \frac{B02}{B04} \\ y=3.896-4.186 \frac{B02}{B03}\)	B02, B03, B04	\(R^{2}\), RMSE	Río4	Landsat-8	[5]
\(y=37661 \cdot B8A^{2}+1845 \cdot B8A <br> y=531.5- \frac{B04}{0.88}\)	B04, B8A	\(R^{2}\), RMSE, MAPE	Estuario5	Pléiades	[6]

Múltiples modelos (lineal, logaritmos, inversa, cuadrática, exponencial, potencial) y plataformas (Sentinel-2, Landsat-5 y Landsat-8) emplean el cociente de bandas B04/B03 [7].

Modelos de estimación a partir de Sentinel-2 y Landsat-8 utilizan regresiones lineales, cuadráticas y logarítmicas empleando B02, B03, B04, B01 (con menos apariciones) y cocientes entre éstas [8].

1.2 Métodos de aprendizaje automático

El aprendizaje automático es un subconjunto de la inteligencia artificial que permite que un sistema aprenda y mejore de forma autónoma, sin necesidad de una programación explícita, a través del análisis de grandes cantidades de datos. El aprendizaje automático permite que los sistemas informáticos se ajusten y mejoren continuamente a medida que acumulan más “experiencias”. Por lo tanto, el rendimiento de estos sistemas puede mejorar si se proporcionan conjuntos de datos más grandes y variados para su procesamiento.

Cuando se entrenan modelos de machine learning, cada conjunto de datos y cada modelo necesitan un conjunto diferente de “hiperparámetros”. Los hiperparámetros son variables de configuración externa que se utilizan para administrar el entrenamiento de modelos de machine learning. Controlan de forma directa la estructura, funciones y rendimiento de los modelos. Los hiperparámetros son los parámetros de un modelo de aprendizaje automático, que no se aprenden durante el entrenamiento, sino que se establecen antes de que comience.

El “ajuste de hiperparámetros” permite modificar el rendimiento del modelo para lograr resultados óptimos. Este proceso es una parte fundamental del machine learning. El ajuste de hiperparámetros puede ser manual o automático. A pesar de que el ajuste manual es lento y tedioso, permite entender mejor cómo afectan al modelo las ponderaciones de los hiperparámetros. El proceso de ajuste de hiperparámetros es iterativo, y debe probar diferentes combinaciones de parámetros y valores.

En el aprendizaje automático es importante utilizar técnicas de “validación cruzada” , de modo que el modelo no se centre únicamente en una única porción de sus datos. La validación cruzada o cross-validation es una técnica utilizada para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba. La idea básica de la validación cruzada es dividir los datos en conjuntos de entrenamiento y validación, y luego entrenar el modelo en el conjunto de entrenamiento y evaluar su rendimiento en el conjunto de validación. Este proceso se repite varias veces, con diferentes subconjuntos de los datos utilizados para el entrenamiento y la validación, y se calcula el rendimiento promedio.

En los procesos de machine learning supervisado se utilizan diversos algoritmos y técnicas de cálculo, generalmente calculados mediante el uso de programas como R o Python.

Dependiendo del tipo de datos que se usen para el entrenamiento, será de modelo de aprendizaje automático que se use. A grandes rasgos, existen tres tipos de modelos que se usan en el aprendizaje automático: aprendizaje supervisado , no supervisado y por refuerzo.

Consultando el trabajo de otros investigadores, se observa que utilizan principalmente el aprendizaje automático supervisado. Este tipo aprendizaje supervisado utiliza un conjunto de entrenamiento para enseñar a los modelos a producir el resultado deseado. Este conjunto de datos de entrenamiento incluye entradas y salidas correctas, que permiten al modelo aprender con el tiempo. El algoritmo mide su precisión a través de la función de pérdida, ajustando hasta que el error se haya minimizado lo suficiente.

Yang Zhe y otros, utilizaron como datos de entrada la reflectancia de superficie y datos de salida la turbidez, utilizaron los modelos SVR (support vector regression), random forest (RF) y eXtreme Gradiente Boostring (XGBoost). Los hiperparámetros de cada modelo se determinaron mediante una búsqueda en cuadrícula de validación cruzada en Scikit-Learn de Python [9].

Ma Yue y otros, utilizaron varios modelos de aprendizaje automático, usaron Python 3.7 tanto para la predicción de la turbidez del agua y la optimización de la los hiperparámetros [2].

Zhao y otros probaron 14 modelos de machine learning en un estanque de peces con un dispositivo de construction propia, de los cuales ETR, Bagging, RFR, and ABR son los que presentaron un mejor desempeño en la estimación de la turbidez. Los algoritmos se implementaron utilizando Python 3.6 y bibliotecas de aprendizaje scikit [10].

Tabla 2: Características principales de algoritmos de aprendizaje automático para la estimación de turbidez.

Modelo de machine learning	Cuerpo de agua	Métricas	Plataforma	Referencia
SVR, ELM ,BP ,CART ,GBT ,RF ,KNN	Lagos	RMSE, \(R^{2}\), MAE	Sentinel-MSI	[2]
eXtreme Gradient Boosting (XGBoost), support vector regression (SVR), random forest (RF)	Lago	RMSE, \(R^{2}\), MAPE	Sentinel-2A/B y Landsat-8/9	[9]
linear regression (LR), ridge regression (RR), least absolute shrinkage and selection operator regression(LASSO), elastic net regression (ENR), k-nearest neighbor regression (KNN), Gaussian process regression (GPR), decision tree regression (DTR), support vector regression (SVR), multilayer perceptron regression (MLP), adaptive boosting regression (ABR), gradient boosting regression (GBR), bootstrap aggregating regression (Bagging), random forest regression (RFR), and extreme tree regression (ETR)	Estanque de peces	MAE, MRSE, MAPE, \(R^{2}\), RE, Acc	Propia	[10]

2 Lectura y procesamiento de datos

Este sitio web está configurado de manera tal que se mantenga actualizado al tener disponibles nuevos datos de entrada, que serán procesados posteriormente. Mostrando siempre los resultados de todo el proceso.

El código a continuación se dedica a la descarga automática de archivos desde el repositorio principal del proyecto.

from descarga_datos import main
main()

2.1 Procesamiento de archivos csv

Para el procesamiento de los datos se utilizará la librería pandas de Python.

En el proyecto tenemos dos archivos .csv que contienen los datos:

-base_de_datos_lab.csv → contiene resultados de laboratorio

-base_de_datos_gis.csv → contiene datos espectrales

import pandas as pd

#Lectura de datos

df1_lab = pd.read_csv(r"D:\GIT\Turbidez\Input\base_de_datos_lab.csv") # DataFrame de datos de laboratorio.
df2_gis = pd.read_csv(r"D:\GIT\Turbidez\Input\base_de_datos_gis_acolite.csv") # DataFrame de datos espectrales

# Combinamos ambos DataFrame para tener los datos en una única base de datos

df_combinado = pd.merge(df1_lab, df2_gis, on=['latitud', 'longitud','fecha'], how='inner')

# Filtramos las filas que contengas los valores de turbidez

df_turbidez = df_combinado[(df_combinado['param'] == 'turb')]

# Eliminamos las columnas que no nos interesan
df_turbidez_banda = df_turbidez.drop(columns=['longitud','latitud','punto','fecha','param'])

df_turbidez_banda.head()

# Cambiamos el nombre de la columna "valor" por el de "turbidez

df_turbidez_banda.rename(columns={'valor': 'turbidez'}, inplace=True)

df_turbidez_banda.head()

# Pivotamos tabla

df_final = df_turbidez_banda.pivot_table(
        index='turbidez',
        columns='banda',
        values='reflect',)

# Creamos tabla final

df_final.to_csv(r"D:\GIT\Turbidez\Output\csv\Turbidez-Bandas.csv", index=True)

La datos útiles a este sitio están ordenanos de la siguiente forma. (Solo se muestran las 5 primeras filas)

	turbidez	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A
0	2.0	0.030264	0.035694	0.052761	0.038619	0.028929	0.017878	0.018247	0.014811	0.008830	0.006162	0.015677
1	3.0	0.018744	0.027450	0.044980	0.028040	0.018781	0.007693	0.007571	0.007623	0.001773	0.000112	0.003305
2	4.0	0.029141	0.036459	0.054631	0.039793	0.029767	0.022395	0.022966	0.019097	0.009430	0.007976	0.020290
3	4.5	0.022306	0.027934	0.045472	0.028727	0.020835	0.014462	0.017479	0.013988	0.000088	0.000361	0.013232
4	5.0	0.029902	0.030043	0.047235	0.033111	0.024060	0.011116	0.012519	0.007006	0.001982	0.000181	0.009074

2.2 Procesamiento de productos satelitales

Para el procesamiento se utilizará la librería rasterio en python.

Tabla de referencia para seleccionar bandas en rasterio

Nombre de banda en archivo .tif	Banda de Sentinel-2	Número para la lectura en rasterio
Banda 1	B01	1
Banda 2	B02	2
Banda 3	B03	3
Banda 4	B04	4
Banda 5	B05	5
Banda 6	B06	6
Banda 7	B07	7
Banda 8	B08	8
Banda 9	B8A	9
Banda 10	B11	10
Banda 11	B12	11

Lectura de imágen raster y visualización del área de estudio, con el objetivo de verificar que el código está funcionando y de que se leen correctamente los archivos.

Se utilizrá el índice NDWI para recortar cada la imágen satelital, y poder trabajar únicamente con los píxeles de agua.

\(NDWI=(B03-B11)/(B03+B11)\)

import numpy as np
import matplotlib.pyplot as plt
from skimage.filters import threshold_otsu

# NDWI
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)

# Calcular umbral automático con Otsu
umbral = threshold_otsu(ndwi[np.isfinite(ndwi)])  # descartamos NaN e inf

# Binarización: agua = 1, no agua = 0
ndwi_bin = np.where(ndwi >= umbral, 1, 0)

# Gráfico máscara de agua
fig, ax = plt.subplots(figsize=(4.5, 4.5))
img = ax.imshow(ndwi_bin, cmap='Greys')
ax.set_title(f"NDWI - Otsu (umbral={umbral:.3f})")
cbar = fig.colorbar(img, ax=ax, fraction=0.046, pad=0.04)
cbar.set_label("Agua (1) / No Agua (0)")
ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
ax.axis("on")
plt.tight_layout()
plt.show()

Durante la verificacion de los recortes de las imágenes satélitales, se observó que con el método Otsu no se puede calcular correctamente el valor umbral agua/no agua , por ende no se puede procesar correctamente algunas imágenes. Para resolver esto se decidió utilizar el valor promedio del NDWI para calcular el valor umbral, el cual sí puede hacer una mejor demarcación de las zona agua/no agua.

Una vez entrenado y validado los modelos, se hará automáticamente la lectura de archivos .tif de entrada, a cada una se realizara el cálculo de umbral de máscara de agua y su procesamiento para generar mapas de turbidez,

2.3 Pruebas de correlación

Para ser más rigurosos, agregaremos ésta etapaa durante el entrenamiento de nuestro modelo lineal.

En esta primera etapa nos interesa estimar la turbidez por medio de modelos lineales sencillos, para encontrar mejores modelos podemos aplicar transformaciones logarítmicas con el objetivo principal de linealizar relaciones no lineales.

Primeramente se probará un modelo sin aplicar ninguna transformaciones en las variables, luego se aplicará logaritmo natural en la turbidez, en las bandas espectrales y en en ambas varaibles, con el fin de evaluar en rendimiento de los modelos y seleccionar los mejores, según sus métricas de desempeño.

Se utilizará el coeficiente de correlación lineal r, para evaluar la relación lineal entre la turdidez y las distintas bandas (con y sin transformación logarítmica).

Calculamos coeficiente de correlacion lineal r entre la turbidez y cada banda con la función .corr de pandas.

Esta medida indica cuanto se relacionan dos variables, puede tomar valores desde -1 a +1:

• +1 correlación lineal perfecta positiva.

• 0 sin correlación.

• -1 correlación lineal perfecta negativa.

Prueba de correlación C1: (sin aplicar logaritmo)

Cálculo de correlación lineal entre la turbidez y las bandas espectrales.

import pandas as pd

Datos= pd.read_csv(r'D:\GIT\Turbidez\Output\csv\Turbidez-Bandas.csv')

bandas = [col for col in Datos.columns if col.startswith('B')]

correlaciones = {}

for banda in bandas:
    r = Datos['turbidez'].corr(Datos[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones1 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones1 = df_correlaciones1.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones1['Banda'] = "turb vs " + df_correlaciones1['Banda'].astype(str)
df_correlaciones1.rename(columns={'Banda': 'Combinación 1'}, inplace=True)

df_correlaciones1.to_csv(r'D:\GIT\Turbidez\Output\correlaciones\C1_turb_vs_banda.csv', index=False)

#df_correlaciones1

Prueba de correlación C2: logaritmo en la turbidez

El código a continuación aplica logaritmo a la columna de turbidez.

import numpy as np

Datos_turb_log = pd.read_csv(r'D:\GIT\Turbidez\Output\csv\Turbidez-Bandas.csv')
Datos_turb_log['turbidez'] = np.log(Datos_turb_log['turbidez'])

#Cambio el nombre la columna "turbidez" luego de aplicar el logaritmo
Datos_turb_log = Datos_turb_log.rename(columns={'turbidez': 'ln_turbidez'})

#Creo archivo csv para usarlo para entrenar modelos
Datos_turb_log.to_csv(r"D:\GIT\Turbidez\Output\csv\ln_turb-bandas.csv", index=False)

Cálculo de correlación lineal entre el logaritmo de la turbidez y las bandas espectrales.

bandas = [col for col in Datos_turb_log.columns if col.startswith('B')]

correlaciones = {}

for banda in bandas:
    r = Datos_turb_log['ln_turbidez'].corr(Datos_turb_log[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones2 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones2 = df_correlaciones2.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones2['Banda'] = "ln_turb vs " + df_correlaciones2['Banda'].astype(str)
df_correlaciones2.rename(columns={'Banda': 'Combinación 2'}, inplace=True)

df_correlaciones2.to_csv(r'D:\GIT\Turbidez\Output\correlaciones\C2_ln_turb_vs_banda.csv', index=False)

Prueba de correlación C3: logaritmo en las bandas

Anteriormente en el proyecto se estaba trabajando con una base de datos de reflectancias procesadas con Sen2Cor. Actualmente trabajamos con una nueva base de datos con las reflectancias corregidas por Acolite.

Acolite y Sen2Cor son dos herramientas de software utilizadas para la corrección atmosférica de imágenes satelitales. Acolite se destaca en ambientes acuáticos eutróficos, mientras que Sen2Cor es un procesador más general para la creación de productos Sentinel-2 de nivel 2A.

IMPORTANTE: Este cambio provocó errores en el código a partir de esta sección en adelante, el error se debe a que hay valores de reflectancia que son cero y negativos en la banda 12, en los datos provenientes de ACOLITE.

En las pruebas de correlación C3 y C4 aplicamos logaritmo a las bandas. Por definción de logaritmo no ponemos tomar valores cero ni negativos. Para salvar este error, se procedió a filtrar filas con esos valores.

Creamos un nuevo DataFrame para aplicarle el logaritmo a las bandas. (Se usará tambien para la prueba de correlación 4)

import pandas as pd

# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Output\csv\Turbidez-Bandas.csv')

# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrados = Datos[(Datos > 0).all(axis=1)]


Datos_filtrados.to_csv(r"D:\GIT\Turbidez\Output\csv\Datos_filtrados-para-usar-ln.csv", index=False)

Aplicamos logaritmo a las bandas

# Este chunk sólo se usa para aplicar logaritmo a datos_filtrados 
import numpy as np

DatosC3 = Datos_filtrados.copy()

col = [col for col in DatosC3.columns if col.startswith('B')]

DatosC3[col] = np.log(DatosC3[col])

#Cambiamos en nombre las columnas, agremamos ln_ a cada columna

DatosC3.columns = ['ln_' + col for col in DatosC3.columns]

DatosC3.rename(columns={'ln_turbidez': 'turbidez'}, inplace=True)

#Creo archivo csv para usarlo para entrenar modelos
DatosC3.to_csv(r"D:\GIT\Turbidez\Output\correlaciones\C3_turb_vs_ln_banda.csv", index=False)

Cálculo de correlación lineal entre la turbidez y el logaritmo de las bandas espectrales.

bandas_ln = [col for col in DatosC3.columns if col.startswith('ln_B')]

correlaciones = {}

for banda in bandas_ln:
    r = DatosC3['turbidez'].corr(DatosC3[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones3 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones3 = df_correlaciones3.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones3['Banda'] = "turb vs " + df_correlaciones3['Banda'].astype(str)
df_correlaciones3.rename(columns={'Banda': 'Combinación 3'}, inplace=True)


df_correlaciones3.to_csv(r'D:\GIT\Turbidez\Output\correlaciones\C3_turb_vs_ln_banda.csv', index=False)

Prueba de correlación C4: logaritmo en turbidez y en bandas

import numpy as np

DatosC4 = np.log(Datos_filtrados)
#Cambiamos en nombre las columnas, agremamos ln_ a cada columna
DatosC4.columns = ['ln_' + col for col in DatosC4.columns]

#Creo archivo csv para usarlo para entrenar modelos
DatosC4.to_csv(r"D:\GIT\Turbidez\Output\csv\ln_turb-ln_banda.csv", index=False)

Cálculo de correlación lineal entre el logaritmo de la turbidez y el logaritmo de las bandas espectrales.

bandas_ln = [col for col in DatosC4.columns if col.startswith('ln_B')]

correlaciones = {}

for banda in bandas_ln:
    r = DatosC4['ln_turbidez'].corr(DatosC4[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones4 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones4 = df_correlaciones4.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones4['Banda'] = "ln_turb vs " + df_correlaciones4['Banda'].astype(str)
df_correlaciones4.rename(columns={'Banda': 'Combinación 4'}, inplace=True)

df_correlaciones4.to_csv(r'D:\GIT\Turbidez\Output\correlaciones\C4_ln_turb_vs_ln_banda.csv', index=False)

Se resume en una tabla las pruebas de correlación realizadas

Tabla resúmen de pruebas de correlación lineal r** **

	Combinación 1	r	Combinación 2	r	Combinación 3	r	Combinación 4	r
0	turb vs B05	0.887299	ln_turb vs B05	0.813367	turb vs ln_B05	0.734057	ln_turb vs ln_B05	0.844506
1	turb vs B06	0.876827	ln_turb vs B04	0.752422	turb vs ln_B06	0.717491	ln_turb vs ln_B04	0.787566
2	turb vs B08	0.869149	ln_turb vs B06	0.714446	turb vs ln_B07	0.706154	ln_turb vs ln_B06	0.770889
3	turb vs B07	0.867741	ln_turb vs B07	0.707941	turb vs ln_B04	0.699797	ln_turb vs ln_B07	0.754235
4	turb vs B04	0.827019	ln_turb vs B08	0.703117	turb vs ln_B08	0.698987	ln_turb vs ln_B08	0.740944
5	turb vs B8A	0.788088	ln_turb vs B8A	0.603662	turb vs ln_B8A	0.641109	ln_turb vs ln_B8A	0.661797
6	turb vs B03	0.673358	ln_turb vs B03	0.514314	turb vs ln_B03	0.598708	ln_turb vs ln_B03	0.533858
7	turb vs B02	0.614315	ln_turb vs B02	0.429268	turb vs ln_B02	0.543950	ln_turb vs ln_B02	0.447893
8	turb vs B01	0.611305	ln_turb vs B01	0.391705	turb vs ln_B01	0.533258	ln_turb vs ln_B01	0.397983
9	turb vs B12	0.429200	ln_turb vs B12	0.242050	turb vs ln_B11	0.352344	ln_turb vs ln_B11	0.250618
10	turb vs B11	0.409447	ln_turb vs B11	0.222510	turb vs ln_B12	0.351269	ln_turb vs ln_B12	0.223253

2.4 Seleción de variables

A partir del análisis con las correlaciones entre turbidez (y su logaritmo) y distintas bandas (y sus logaritmos). Se puede observar que en las combinaciones 1 y 4 las correlaciones son mas altas respecto a C2 y C3. Por que que se seguirá el análisis con estas combinaciones:

• Combinación 1: turb vs bandas

• Combinación 4: ln_turb vs ln_bandas

Se entrenarán modelos con ambas combinaciones por separado y se evalurá su desempeño.

Todos los modelos serán sometidos a bootstrapping con enfoque out-of-bag (OOB):

→ Donde se tomará el 75% de los datos con reemplazo para entrenamiento; → Luego se usa el 25% restante (no muestreados) como conjunto de prueba (test).

Proponemos avanzar hacia un modelo multivariable, utilizando el AIC [11] como criterio de seleccion de variables, este criterio es muy conveniente para encontrar el modelo de mejor balance entre ajuste y complejidad.

¿Qué hace el AIC?

→ Penaliza modelos con demasiadas variables para evitar sobreajuste.

→ Se busca minimizar el AIC: cuanto más bajo, mejor el modelo.

→ No requiere que las variables tengan alta correlación individual, pero sí que mejoren el modelo conjunto.

AIC=n⋅ln(RMSE²)+2k

2.4.1 Combinación 1: turb vs bandas

	variables	num_variables	R²_train (bootstrap)	RMSE_train (bootstrap)	R²_test	R²-ajustado	RMSE_test	AIC
0	B05	1	0.817394	35.382538	0.301074	0.257391	33.316077	130.217443
1	B05, B06	2	0.818716	35.254222	0.326811	0.237053	32.696901	131.542091
2	B05, B06, B08	3	0.826294	34.509493	0.270486	0.114162	34.037282	134.988433
3	B05, B06, B08, B07	4	0.842965	32.811820	0.412432	0.231642	30.546900	133.093476
4	B05, B06, B08, B07, B04	5	0.928830	22.089167	0.510597	0.306679	27.878628	131.802973
5	B05, B06, B08, B07, B04, B8A	6	0.942895	19.786427	0.756598	0.623833	19.660758	121.230488
6	B05, B06, B08, B07, B04, B8A, B03	7	0.947098	19.044359	0.725304	0.533016	20.886450	125.407622
7	B05, B06, B08, B07, B04, B8A, B03, B02	8	0.971206	14.050324	0.861429	0.738255	14.834547	115.090513
8	B05, B06, B08, B07, B04, B8A, B03, B02, B01	9	0.976194	12.775463	0.905189	0.798526	12.270664	110.259610
9	B05, B06, B08, B07, B04, B8A, B03, B02, B01, B12	10	0.979567	11.835854	0.923241	0.813586	11.040831	108.457610
10	B05, B06, B08, B07, B04, B8A, B03, B02, B01, B...	11	0.980154	11.664591	0.916229	0.762648	11.534155	112.031255

Gráfico de AIC

¿Por qué el AIC baja al agregar más variables, si se supone que penaliza la complejidad?

AIC=n⋅ln(RMSE²)+2k

donde:

n = número de observaciones k = número de parámetros (incluye el intercepto) ln(RMSE²) representa el ajuste del modelo

AIC no busca modelos simples en sí, sino la mejor compensación entre ajuste y complejidad.

A pesar de agregar más variables, el error (RMSE) está disminuyendo mucho más rápido de lo que penaliza el AIC con +2k

2.4.2 Combinación 4: ln_turb vs ln_bandas

	variables	num_variables	R²_train (bootstrap)	RMSE_train (bootstrap)	R²_test	R²-ajustado	RMSE_test	AIC
0	ln_B05	1	0.683069	0.731393	0.777470	0.752744	0.593922	-7.462168
1	ln_B05, ln_B04	2	0.814241	0.559944	0.768898	0.711123	0.605252	-5.046416
2	ln_B05, ln_B04, ln_B06	3	0.911963	0.385479	0.755483	0.650690	0.622572	-2.425724
3	ln_B05, ln_B04, ln_B06, ln_B07	4	0.909252	0.391369	0.732721	0.554535	0.650905	0.553385
4	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08	5	0.908096	0.393855	0.728706	0.457411	0.655776	2.717401
5	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A	6	0.943807	0.307971	0.561432	-0.096420	0.833784	10.000812
6	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A...	7	0.944500	0.306065	0.581122	-0.396261	0.814852	11.495533
7	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A...	8	0.946288	0.301097	0.645448	-0.772761	0.749679	11.661562
8	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A...	9	0.945395	0.303589	0.731296	-1.687039	0.652637	10.611861
9	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A...	10	0.948380	0.295175	0.727397	NaN	0.657355	12.770329
10	ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A...	11	0.950492	0.289072	0.705746	3.942544	0.682962	15.611041

Gráfico de AIC

3 Entrenamiento de modelo lineal

A continuación se presentan los modelos lineales candidatos para estimar la turbidez, cada uno con sus respectivas métricas de desempeño.

3.1 Modelos lineales

3.1.1 Modelo lineal 1

Utiliza B05, para estimar tubidez

\(\displaystyle turbidez = 1290.9721 \cdot B05 + -39.8799\)

\(\displaystyle turbidez = 1290.9721 \cdot B05 + -39.8799\)

3.1.2 Modelo lineal 2

Utiliza 6 bandas para estimar la turbidez: B04, B05, B06, B07, B08 y B8A

\(\displaystyle turbidez = 434.6144 \cdot B05 + 5574.6016 \cdot B06 + 1796.9343 \cdot B08 + -1571.8740 \cdot B07 + -1151.2411 \cdot B04 + -3930.0577 \cdot B8A + 2.1117\)

3.1.3 Modelo lineal 3

Utiliza el logaritmo de la turbidez y en la banda B05

Escala logarítmica

\(\displaystyle \ln(\text{turbidez}) = 1.5682 \cdot ln(B05) + 7.5989\)

Escala real

3.1.4 Modelo lineal 4

Utiliza el logaritmo de la turbidez y en las bandas: B04, B05 y B06

Escala Logarítmica

\(\displaystyle \ln(\text{turbidez}) = 8.3951 \cdot ln(B05) + -5.4853 \cdot ln(B04) + -1.8377 \cdot ln(B06) + 6.1731\)

\(\displaystyle \text{turbidez} = 479.6568 \cdot B05^{8.3951} \cdot B04^{-5.4853} \cdot B06^{-1.8377}\)

Escala Real

3.1.5 Modelo lineal elegido

Luego del anális de las distantas pruebas de correlación con diferenctes bandas y combinaciones de ellas, haciendo un balance entre complejidad del modelo y métricas de desempeño, se opta por elegir el Modelo lineal 1, que utiliza la B05 como variable predictora de la turbidez y será comparado con algún método de aprendizaje automático próximamente.

Si bien el Modelo lineal 2 (que usa 6 bandas) tiene mejores métricas de desempeño, la incorporacion de tantas bandas lo hace muy complejo. Los modelos que usan transformación logaritmica tienen un peor desempeño, respecto a los modelos sin transformaciones.

A continuación se observa la ecuación del modelo, métricas de desempeño y mapas de turbidez:

Modelo con B05 y sus métricas

\(\displaystyle turbidez = 1292.2227 \cdot B05 + -39.9463\)

\(\displaystyle turbidez = 1292.2227 \cdot B05 + -39.9463\)

Mapas de turbidez modelo lineal elegido

Se procesaron 7 archivos .tif

4 Entrenamiento de modelo de aprendizaje automático

4.1 Modelo XGBoost

Se seleccionó XGBoost, como modelo de aprendizaje automático para la estimación de la turbidez. Se importará el modelo desde sklearn, para su entrenamiento y conocer sus hiperparámetros. En esta etapa se utilizará XGBoost sin optimizacion de hiperparámetros, en una etapa posterior se elegirá un método de optimizacion y se evalurá nuevamente las métricas de desempeño.

Entrenamiento

from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np
import os

# Leer datos

df = pd.read_csv(r"D:\GIT\Turbidez\Output\csv\Turbidez-Bandas.csv")

X = df[["B01","B02","B03","B04","B05","B06","B07","B08","B11","B12","B8A"]]
y = df["turbidez"]

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)

# Modelo XGBoost SIN optimizar

model = XGBRegressor(
n_estimators=500,
learning_rate=0.05,
max_depth=6,
subsample=0.8,
colsample_bytree=0.8,
random_state=42)

model.fit(X_train, y_train)

# Guardar modelo

ruta = r"D:\GIT\Turbidez\Output\Modelos"
os.makedirs(ruta, exist_ok=True)
model.save_model(os.path.join(ruta, "modelo_xgboost_sin_optimizar.json"))

Métricas de desempeño

4.2 Hiperparámetros utilizados

• n_estimators=500 → número de arboles que se entrenan secuencialmente

• learning_rate=0.05 → indica que tanta aprende cada arbol

Es el “paso” con que el modelo corrige los errores. -Valores bajos (0.01–0.1) → aprendizaje más lento pero más estable. -Valores altos (>0.3) → puede sobreajustar rápido.

• max_depth=6 → Profundidad máxima de cada árbol.

Controla la complejidad de las interacciones entre las variables (bandas). Más profundo → aprende patrones complejos, pero también ruido.

• subsample=0.8 → Fracción de los datos usados para entrenar cada árbol.

Introduce aleatoriedad. Ayuda a reducir el sobreajuste, 0.8 significa que cada árbol usa el 80% de las filas del dataset.

• colsample_bytree=0.8 → Fracción de las variables (columnas) usadas para cada árbol. Igual que subsample, pero sobre las bandas.

Esto fuerza a los árboles a usar combinaciones distintas de bandas, lo que mejora la generalización.

0.8 = buena práctica para evitar que el modelo dependa demasiado de una sola banda.

Los hiperparámetros controlan el comportamiento del modelo, pero no se aprenden automáticamente del conjunto de datos. Por eso, necesitamos un método que busque las mejores combinaciones posibles para obtener el máximo rendimiento sin sobreajustar.

4.3 Optimización de hiperparámetros

Existen muchos métodos de optimizacion de hiperparámetros

1. Grid Search (búsqueda en rejilla)

Prueba todas las combinaciones posibles de los valores que uno defina.

Por ejemplo, si se prueba 3 valores de max_depth y 3 de learning_rate, va a entrenar 3 × 3 = 9 modelos.

Usa validación cruzada (cross-validation) para evaluar cada combinación.

Ventajas

• Simple, exhaustivo y fácil de entender.

• Ideal cuando el número de combinaciones es pequeño.

Desventajas

• Muy lento si probás muchos parámetros o valores.

• No aprende de los resultados anteriores (prueba todo a ciegas).

2. Random Search (búsqueda aleatoria)

En lugar de probar todas las combinaciones, elige al azar un número fijo de combinaciones.

Ejemplo: podrías probar 30 combinaciones elegidas aleatoriamente de un rango mucho más grande.

Ventajas

• Mucho más rápido que Grid Search.

• Puede explorar un espacio grande de parámetros.

• Buen punto medio entre precisión y costo computacional.

Desventajas

• Resultados algo aleatorios (aunque reproducibles con random_state).

• No garantiza probar la “mejor” combinación posible.

3. Búsqueda bayesiana inteligente (Optuna / Hyperopt / scikit-optimize)

Usa inteligencia adaptativa: aprende de las pruebas anteriores y dirige las siguientes pruebas hacia las regiones más prometedoras.

Es una búsqueda inteligente, no aleatoria ni exhaustiva.

Ventajas

• Mucho más eficiente que Grid o Random Search.

• Encuentra buenos resultados con menos intentos.

• Ideal para modelos complejos como XGBoost.

Desventajas

• Requiere instalar librerías externas (como optuna).

• Más técnica, aunque muy potente.

En este trabajo se optó por utilar Optuna para la optimizacion de los hiperparámetros

4.4 Modelo XGBoost optimizado con Optuna (Búsqueda bayesiana inteligente)

Tabla con hiperparámetros optimizados

	max_depth	learning_rate	subsample	colsample_bytree	n_estimators	R2_promedio_CV
0	9	0.042543	0.751347	0.892266	649	0.924281

Métricas de desempeño

Mapas XGBoost Optuna

Se procesaron 7 archivos .tif

4.5 Relevancia de Bandas espectrales en el modelo XGBOOST

La importancia de las bandas indica cuánto contribuye cada banda espectral (o combinación de ellas) al rendimiento del modelo. Esto permite saber cuáles longitudes de onda son más informativas para estimar la variable objetivo.

Interpretación:

Cuanto mayor sea la importancia, más influye esa banda en las divisiones del modelo.

No necesariamente significa correlación directa (puede ser no lineal).

Si una banda tiene muy baja importancia (≈0), puede eliminarse en futuras versiones del modelo para simplificarlo.

Cuando entrenás un modelo como XGBoost, el algoritmo aprende qué variables (bandas espectrales) son más útiles para predecir la turbidez.

• La importancia de las bandas indica cuánto aportó cada banda al modelo para mejorar sus predicciones.

En XGBoost, el atributo model.feature_importances_ mide la importancia promedio de cada variable dentro de todos los árboles del modelo.

XGBoost da la importancia relativa (no absoluta ni porcentual) basada en “weight”. Los valores se normalizan para que la suma de todos sea 1.0.

Tipo de importancia (importance_type)	Qué mide exactamente	Interpretación práctica
"weight" (por defecto)	Número de veces que la variable fue usada para dividir nodos	Bandas usadas más veces para tomar decisiones
"gain"	Promedio de la mejora en la predicción al usar esa banda	Cuánto “valor” aporta cada vez que se usa
"cover"	Promedio de la cantidad de datos afectados por esas divisiones	Qué tan general es su influencia
"total_gain" o "total_cover"	Suma total de las métricas anteriores	Importancia acumulada

Importancia de cada banda dentro del modelo XGBoost

	Banda	Importancia
0	B05	0.594959
1	B06	0.146968
2	B01	0.055881
3	B02	0.053704
4	B08	0.046594
5	B8A	0.046546
6	B07	0.037112
7	B03	0.009663
8	B04	0.003476
9	B11	0.002844
10	B12	0.002254

El análisis de importancia de características del modelo XGBoost optimizado muestra una dominancia muy marcada de la banda B05, que concentra cerca del 60% de la importancia total. Esto indica que, para el modelo, la información contenida en B05 es el factor más determinante para estimar la turbidez.

En segundo lugar aparece la banda B06, con una importancia mucho menor pero aún relevante aproximadamente 15%. Después de estas dos, el resto de las bandas (B01, B02, B08, B8A, B07) aportan importancias moderadas y similares entre sí, lo que indica que contribuyen al modelo pero de forma secundaria.

Por otro lado, algunas bandas como B03, B04, B11 y B12 tienen una importancia prácticamente nula, lo que sugiere que su aporte al modelo es mínimo. Esto podría deberse a que estas bandas no contienen información útil para diferenciar niveles de turbidez o dicha información está fuertemente correlacionada con otras bandas más relevantes.

El modelo depende fuertemente de un conjunto reducido de bandas, especialmente B05.

La alta concentración de importancia en pocas bandas sugiere que la señal espectral relacionada con la turbidez está localizada en rangos específicos del espectro.

Las bandas con importancia mínima podrían ser candidatas para eliminación en futuros experimentos, lo que implicaría: -Modelos más ligeros -Menor costo computacional -Potencialmente mejorar generalización al reducir ruido

5 Conclusión

El análisis realizado demuestra que la estimación de turbidez a partir de bandas espectrales Sentinel-2 presenta un buen potencial predictivo, pero la capacidad del modelo depende fuertemente del tipo de algoritmo utilizado.

El modelo lineal seleccionado inicialmente permitió explorar la relación entre cada banda y la turbidez, mostrando que ciertas longitudes de onda —especialmente B05, B06 y B08— concentran la mayor parte de la información relevante. Sin embargo, el comportamiento real del fenómeno es notablemente no lineal, con interacciones entre bandas que un modelo lineal no logra capturar completamente. Esto se reflejó en un desempeño aceptable pero limitado, donde el modelo tendió a subestimar valores altos de turbidez y a producir una mayor dispersión en las predicciones.

En contraste, el XGBoost optimizado mediante Optuna mostró un salto significativo en capacidad predictiva. El proceso de optimización permitió ajustar hiperparámetros clave, logrando un modelo más robusto, con mejor generalización y con un ajuste más preciso tanto para los datos de entrenamiento como para los de testeo. Los indicadores R² y RMSE mejoraron sustancialmente respecto del modelo lineal, y esto se evidencia visualmente en el gráfico estimado vs real, donde la dispersión es mucho menor y la tendencia se ajusta mejor a la diagonal ideal.

Además, el análisis de importancia de características del XGBoost confirma y profundiza lo detectado con el modelo lineal: las bandas del espectro rojo y NIR (especialmente B05) son determinantes para explicar la turbidez. Sin embargo, el modelo no lineal también asigna pequeñas contribuciones a otras bandas que, aunque individualmente débiles, aportan información útil cuando interactúan entre sí, lo que el modelo lineal no puede aprovechar.

El modelo lineal es simple, interpretable y útil para comprender relaciones básicas entre bandas y turbidez, pero su precisión es limitada.

El XGBoost optimizado ofrece una mejora clara en desempeño, captura relaciones no lineales y reduce significativamente el error, mostrando ser el enfoque más adecuado para el problema.

La combinación de optimización automática (Optuna) y un algoritmo avanzado de ensamble permite obtener un modelo robusto y reproducible, adaptable a distintos conjuntos de datos y con potencial para implementarse operativamente.

En conclusión, el modelo XGBoost optimizado se posiciona como la mejor alternativa para estimar turbidez a partir de productos satelitales Sentinel-2, superando ampliamente al modelo lineal tanto en precisión como en capacidad de generalización.

Referencias

[1]

J. cc Delegido et al., «Turbidity and secchi disc depth with sentinel-2 in different trophic status reservoirs at the comunidad valenciana», Revista de Teledeteccion, vol. 2019, pp. 15-24, 2019, doi: 10.4995/raet.2019.12603.

[2]

Y. Ma et al., «Remote sensing of turbidity for lakes in Northeast China using sentinel-2 images with machine learning algorithms», IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 14, pp. 9132-9146, 2021, doi: 10.1109/JSTARS.2021.3109292.

[3]

A. K. M. A. Hossain, C. Mathias, y R. Blanton, «Remote sensing of turbidity in the tennessee river using landsat 8 satellite», Remote Sensing, vol. 13, sep. 2021, doi: 10.3390/rs13183785.

[4]

M. Najafzadeh y S. Basirian, «Evaluation of River Water Quality Index Using Remote Sensing and Artificial Intelligence Models», Remote Sensing, vol. 15, may 2023, doi: 10.3390/rs15092359.

[5]

M. Allam, M. Y. A. Khan, y Q. Meng, «Retrieval of turbidity on a spatio-temporal scale using landsat 8 SR: A case study of the Ramganga River in the Ganges Basin, India», Applied Sciences (Switzerland), vol. 10, jun. 2020, doi: 10.3390/app10113702.

[6]

Y. Luo, D. Doxaran, y Q. Vanhellemont, «Retrieval and validation of water turbidity at metre-scale using Pleiades satellite data: A case study in the Gironde Estuary», Remote Sensing, vol. 12, mar. 2020, doi: 10.3390/rs12060946.

[7]

M. Shen, S. Wang, Y. Li, M. Tang, y Y. Ma, «Pattern of turbidity change in the middle reaches of the yarlung zangbo river, southern tibetan plateau, from 2007 to 2017», Remote Sensing, vol. 13, pp. 1-25, ene. 2021, doi: 10.3390/rs13020182.

[8]

Y. O. Ouma, K. Noor, y K. Herbert, «Modelling Reservoir Chlorophyll- a, TSS, and Turbidity Using Sentinel-2A MSI and Landsat-8 OLI Satellite Sensors with Empirical Multivariate Regression», Journal of Sensors, vol. 2020, 2020, doi: 10.1155/2020/8858408.

[9]

Z. Yang et al., «Combined Retrievals of Turbidity from Sentinel-2A/B and Landsat-8/9 in the Taihu Lake through Machine Learning», Remote Sensing, vol. 15, sep. 2023, doi: 10.3390/rs15174333.

[10]

Y. Zhao et al., «Retrieval of Water Quality Parameters Based on Near-Surface Remote Sensing and Machine Learning Algorithm», Remote Sensing, vol. 14, nov. 2022, doi: 10.3390/rs14215305.

[11]

S. Chatterjee y A. S. Hadi, Regression analysis by example. John Wiley & Sons, 2015.

Notas

1. 0,83 - 112,26 NTU.

2. 2,3 - 107,02 NTU.

3. IOA = index of agreement
   SI = scatter index.

4. 20,6 - 112 NTU
   2,3 - 15,4 NTU.

5. MAPE = Mean Absolute Percentage Error
   0 - 1300 NTU
   0 - 80 NTU.