Sitio de estudio

La presente investigación se realizó utilizando datos de estaciones meteorológicas automatizadas (EMAs) ubicadas en el estado de Campeche, México (Figura 1). El clima predominante es cálido subhúmedo, que se presenta en el 92 % de su territorio y el 7.75 % presenta clima cálido húmedo localizado en la parte este del estado y en la parte norte un porcentaje del 0.05 % con clima semi-seco. La temperatura más alta es mayor a 30 °C y la mínima de 18 °C. La temperatura media anual es de 26 a 27 °C. Las lluvias son de abundantes a muy abundantes durante el verano. La precipitación total anual varía entre 1 200 y 2 000 mm, y en la región norte, de clima semi-seco, es alrededor de 800 mm anuales (^{INEGI, 2017}).

Figura 1 Ubicación de las estaciones meteorológicas en el estado de Campeche, México, usadas en este estudio. 

Las bases de datos cada 10 minutos se obtuvieron de las estaciones meteorológicas automatizadas de la Comisión Nacional del Agua (Conagua) de México. La Tabla 1 muestra la información geográfica de las estaciones meteorológicas utilizadas en este estudio, así como los periodos de tiempo de registro de cada estación meteorológica.

Manejo de datos faltantes y calidad

Las bases de datos fueron procesadas cada 10 minutos, detectando mediante algoritmos implementados en el software Microsoft Excel® series de tiempo faltantes, las cuales posteriormente fueron rellenadas utilizando la técnica de interpolación llamada interpolación polinómica de Hermite (PCHIP). Para una descripción más detallada, consultar ^{Salazar, Ureña y Gallego (2010)}, y ^{Torrente-Cantó (2018)} . Una vez que las series de datos fueron completadas, se construyeron bases de datos diarias. Asimismo, los datos se analizaron para identificar valores atípicos, donde aquellos valores por encima de tres desviaciones estándar de la media se marcaron como atípicos. Se analizaron los datos señalados como atípicos: si un valor atípico extremo inferior se encontraba asociado con un evento de lluvia no se eliminaba; en caso contrario, se eliminaron; esto, con el objetivo de tener modelos funcionales incluso en época lluviosa. En el caso de los valores atípicos extremos superiores fueron eliminados. En ambos casos se utilizó la técnica de interpolación PCHIP para su relleno. La Figura 2 muestra el gráfico de caja y bigotes de las estaciones meteorológicas implicadas en los modelos de este estudio. El bigote representa el mínimo y el máximo de las variables. Respecto a la temperatura máxima, la media varía entre 33 y 34 °C, con valores máximos entre 42 y 43 °C (excepto la estación de Cd. del Carmen); el valor mínimo de la temperatura máxima varía entre 22 y 25 °C. La temperatura mínima tiene una media entre 20 y 24 °C, con valores máximos entre 25 y 28 °C, y mínimos entre 12 y 17 °C. La Figura 2 también muestra los valores atípicos, generalmente asociados con la temperatura mínima.

Figura 2 Gráficos de caja y bigotes de variable temperatura en las estaciones meteorológicas analizadas. 

Ecuación de la FAO 56 PM (ETo-FAO56PM)

La ecuación FAO56PM es el modelo estándar usado para estimar con precisión la ETo, propuesto por la Organización de las Naciones Unidas para la Agricultura y la Alimentación (FAO, por sus siglas en inglés). Incorpora aspectos termodinámicos y aerodinámicos, y considera muchos parámetros meteorológicos relacionados con el proceso de evapotranspiración, como la radiación neta, la temperatura del aire, el déficit de presión de vapor, la velocidad del viento. Ha demostrado ser un método relativamente preciso bajo diferentes condiciones o regiones (^{Allen et al., 1998}). Estos aspectos se han incorporado en la siguiente ecuación:

donde R_n = radiación neta en la superficie (MJ m^-2 día^-1); G = flujo del calor de suelo (MJ m^-2 día^-1); T_med = temperatura media del aire a 2 m de altura (°C); u₂ = velocidad del viento a 2 m de altura (m s^-1); e_s = presión de vapor de saturación (kPa); e_a = presión real de vapor (kPa); ∆ = pendiente de la curva de presión de vapor (kPa °C^-1); γ = constante psicrométrica (kPa °C^-1).

En este estudio, el método ETo-FAO56PM se usó para evaluar los métodos convencionales y de inteligencia artificial.

Ecuación de Hargreaves y Samani

El modelo de HS se considera como un modelo alternativo para estimar la ETo cuando solo los registros de temperatura están disponibles en el lugar de estudio; es uno de los métodos que ha sido utilizado consecutivamente por su simple implementación y la precisión en sus resultados (^{Gong et al., 2016}; ^{Shiri, 2017}). La Ecuación (2) del modelo de Hargreaves y Samani está estructurada de la siguiente manera:

donde ETo = evapotranspiración de referencia (mm día^-1); K_HG = es un coeficiente empírico, que inicialmente fue establecido en 0.0023, pero se ha recalibrado acorde con el lugar empleado; T_med = temperatura media; T_máx = temperatura máxima; T_mín = temperatura mínima; R_a = radiación solar extraterrestre. R_a se calculó en función del día del año, latitud del sitio y ángulo solar de acuerdo con la ecuación propuesta por ^{Allen et al. (1998)}.

Ecuación de Camargo

El modelo de Camargo es una modificación a la ecuación de Thornthwaite (TH). Es un modelo basado en la variable climática de temperatura. Camargo sustituyó el valor de la temperatura media de la ecuación de Thornthwaite por la temperatura media efectiva (T_ef) (Camargo, Marin, Sentelhas, & Picini, 1999). La Ecuación (3) del modelo de Camargo está estructurada de la siguiente manera:

donde ETo = evapotranspiración de referencia (mm día^-1); K_CA1 y K_CA2 = coeficientes empíricos, donde sus valores originales son 16 y 0.36, respectivamente, y deben ser calibrados acorde con el lugar de empleo; I = índice de calor anual; a = exponente empírico en función de I; N = tiempo máximo de insolación en horas; (K_CA2 * (3T_máx - T_mín) = temperatura efectiva, reemplazando a la temperatura media en la ecuación de Thornthwaite. El valor de I se define como la suma de 12 valores de índices de calor mensual, como se muestra en la siguiente ecuación:

donde T_medj = temperatura media mensual (°C).

y:

El valor de a varía de 0 a 4.25; mientras que el índice de calor anual I varía de 0 a 160.

Ajuste de parámetros en los métodos convencionales

Los métodos convencionales basados en temperatura deben ajustarse a las condiciones locales antes de ser utilizados (^{Almorox et al., 2018}) para obtener buenas estimaciones de ETo. Por tal motivo, los coeficientes originales de las ecuaciones se calibraron utilizando técnicas de regresión no lineal mediante al algoritmo Levenberg-Marquardt.

Máquinas de soporte vectorial (SVM)

La técnica máquinas de soporte vectorial (SVM) fue desarrollada por ^{Vapnik (2000)} y es uno de los enfoques basado en el aprendizaje automático. Es una técnica de aprendizaje supervisado, robusta, para resolver problemas de clasificación y regresión aplicados a grandes conjuntos de datos complejos con ruido; selecciona un hiperplano único de separación de cada clase y la idea básica es mapear los datos x en un espacio de características de alta dimensión a través de un mapeo no lineal y hacer una regresión lineal en este espacio. Durante el entrenamiento solo se consideran los ejemplos que se encuentran en el margen de separación, llamados vectores de soporte. Se aplican con éxito en problemas de regresión, generalmente denominados SVR (regresión de vectores de soporte), utilizando SVM para un conjunto de datos {(X_i , Y_i )} N/i = 1, donde X_i es el vector de entrada; Y_i , el valor de salida, y N es el número total de conjuntos de datos mediante el mapeo de X en un espacio característico a través de una función no lineal φ(x) para después encontrar una función de regresión (^{Fan et al., 2018a}; ^{Mehdizadeh, Behmanesh, & Khalili, 2017}; ^{Quej, Almorox, Arnaldo, & Saito, 2017}; ^{Topi & Vanita, 2017}; ^{Wen et al., 2015}):

donde φ (x) es la función de mapeo no lineal; ω es un vector de peso, y b es un valor de sesgo; son los parámetros de la función de regresión, los cuales pueden ser calculados minimizando la siguiente función de riesgo regularizado:

donde el término 12 ω² mejora la generalización del modelo SVM, normalizando el grado de complejidad del modelo. C es un parámetro de compensación positiva que determina el grado de error en el problema de optimización elegido por el usuario; (ε) es la función de pérdida de Vapnik (tamaño del tubo del modelo de SVM) y está definida como:

Es decir, si la diferencia entre los valores predichos y los medidos es menor que ε, entonces la pérdida es igual a 0. Si los valores predichos están dentro del tubo, el error de pérdida es igual a 0. Para el resto de los puntos predichos encontrados fuera del tubo, la pérdida es igual a la diferencia entre el valor predicho y el radio ε del tubo. Para la detección de valores atípicos, las variables de holgura ξ y ξ ⃰ miden de arriba y abajo en el tubo de ε.

Debido a que ambas variables adquieren valores positivos, se tiene que minimizar el riesgo con la siguiente ecuación:

Sujeto a yi - ωϕxi - bi ≤ ε + ξi ωϕxi+ bi - yi ≤ε + ξi*ξ, ξi*≥0 [/p]

donde la C∑i=1n(ξi + ξi*) controlan los grados de riesgo empírico.

El modelo SVM estima la regresión en función de una serie de funciones Kernel, que convierten los datos de entrada originales de dimensiones inferiores a un espacio de características de mayor dimensión de una manera implícita. Entre los Kernel más utilizados se halla el SVM polinomial (SVM-Poly) y la función de base radial SVM (SVM-RBF), cuyos parámetros del Kernel deberán ajustarse previamente mediante un algoritmo. Por ejemplo, los parámetros óptimos del Kernel y del modelo SVM generalmente se obtienen utilizando el método de búsqueda de cuadrícula (^{Mehdizadeh et al., 2017}):

Implementación del modelo SVM en la estimación de la ETo

En el presente estudio, el modelo SVM para estimar la ETo se construyó usando el software R (^{RDevelopment, 2009}). Como variables de entrada se usaron los datos meteorológicos de T_máx, T_mín y R_a , y como variable objetivo se utilizaron los valores de ETo-FAO56PM (Ecuación (1)). Para el entrenamiento y validación del modelo SVM, se utilizó el software R en conjunto con el paquete LIBSVM 3.1 (^{Chang, Lin, & Tieleman, 2013}). Para la redimensión de los datos se utilizó la función Kernel de Base Radial (RBF) (Ecuación (10)). Con la finalidad de evitar el sobre ajuste y aumentar el desempeño del modelo SVM para estimar la ETo, los parámetros ε, C y ɣ de la SVM, y de la función Kernel de base radial se optimaron mediante el algoritmo genético (GA), utilizando validación cruzada (VC = 5 fólders) (^{Quej et al., 2017}; ^{Shrestha & Shukla, 2015}), y variando los parámetros ε = 0.002 a ε = 2, C = 0.0001 a C = 10, y ɣ = 0.0001 a ɣ = 2. El GA se implementó utilizando el software R en conjunto con la librería e1079 y Caret. Se utilizó el 60 % de los datos durante la etapa de entrenamiento y el 40 % en la etapa de validación.

Los parámetros optimizados por GA utilizados en el entrenamiento de la SVM se muestran en la Tabla 2.

Programación de expresión genética (GEP)

La programación de la expresión génica (GEP) fue presentada por ^{Ferreira (2001)} . Es una rama de los algoritmos evolutivos que tiene la capacidad de modelar los procesos dinámicos y no lineales. Es un algoritmo que pertenece a la familia de los algoritmos genéticos (GA) y programación genética (GP) tradicionales. Puede emular la evolución biológica basada en la programación por computadora para resolver un problema definido por el usuario.

Los GEP se consideran un híbrido entre los GA y GP. Utilizan programación genética para la solución del problema en forma de árbol, donde existen dos tipos de nodos:

En GEP, los individuos se codifican primero como cadenas lineales de longitud fija como en GA. Luego, se expresan como entidades no lineales de diferentes tamaños y formas, como el GP (^{Ferreira, 2001}). Además, un conjunto de terminales (coeficientes y predictores), funciones y operadores matemáticos se utilizan en el GEP para estimar la variable dependiente (^{Mehdizadeh et al., 2017}), creando funciones de manera aleatoria y seleccionando aquellas que presenten un mejor ajuste a los resultados experimentales, permitiendo la generación de algoritmos y expresiones matemáticas de manera automática para la solución de problemas (^{Mattar, 2018}; ^{Shiri, 2017}).

Implementación del modelo GEP en la estimación de la ETo

En este estudio, la implementación de la técnica GEP se llevó a cabo utilizando el software GenexprooTols v. 5.0. Las variables de entrada son los valores de los datos meteorológicos de T_máx, T_mín, R_a y valores de ETo-FAO56PM como variable objetivo. Los operadores aritméticos y funciones matemáticas implementadas dentro del programa fueron {+,-, ×, ÷, x, x3,  x2, x3, Inx, eᵡ,sin⁡x,cos⁡x, Arctanx}, recomendados para estudios hidrológicos (^{Mattar, 2018}; ^{Mehdizadeh et al., 2017}; ^{Shiri, 2017}). El 70 % de los datos se usaron para la etapa de entrenamiento y el 30 % para la validación, utilizando validación cruzada (VC = 5 fólders) para evitar el sobreajuste. Los parámetros GEP utilizados en el presente estudio se muestran en la Tabla 3 (^{Shiri et al., 2014}).

XGBoost

Es uno de los algoritmos más importantes y potentes del “Machine Learning” (aprendizaje automático) creado por ^{Chen y Guestrin (2016)}, utilizado para el análisis de problemas de regresión y clasificación estadística, el cual produce un modelo de predicción complejo a partir del ensamblaje de árboles de decisión (modelos simples) en un contexto de aprendizaje supervisado.

El modelo se basa en la teoría del aumento de gradiente, por lo que las predicciones de varios aprendices "débiles" (modelos cuyas predicciones son ligeramente mejores que las suposiciones aleatorias) se combinan para desarrollar un aprendiz "fuerte". Estos aprendices "débiles" se combinan siguiendo una estrategia de aprendizaje gradual, evitando un sobreajuste y optimizando los recursos de cómputo. Esto se obtiene simplificando todas las funciones que permitan combinar términos predictivos y de regularización, pero que, a su vez, mantenga una velocidad computacional óptima durante todo el procesamiento. Al comienzo del proceso de calibración, un aprendiz "débil" se ajusta a todo el espacio de datos, y luego un segundo aprendiz se ajusta a los residuos de la primera. Este proceso de ajuste de un modelo a los residuos del anterior continúa hasta que se alcanza algún criterio de detención (minimización de la raíz cuadrada media del error). El resultado es un tipo de media ponderada de las predicciones individuales de cada alumno débil. Tradicionalmente, los árboles de regresión se seleccionan como aprendices "débiles" (^{Fan et al., 2018a}). Bajo este contexto el modelo XGBoost se basa en la siguiente función objetivo: pérdida + regularizador:

donde l es el término predictivo y Ω el término de regularización. La función de pérdida para el término predictivo puede ser especificada por el usuario. El término de regularización se obtiene con una expresión analítica basada en el número de hojas del árbol y las puntuaciones de cada hoja. El punto clave del proceso de calibración de XGBoost es que ambos términos se reordenan en última instancia en la siguiente expresión:

donde G y H se obtienen de la expansión de las series de Taylor de la función de pérdida; λ es el parámetro de regularización, y T es el número de hojas en un árbol. Esta expresión analítica de la función objetivo permite un rápido escaneo de izquierda a derecha de las posibles divisiones del árbol, pero siempre teniendo en cuenta la complejidad.

XGBoost tiene una amplia gama de parámetros de ajuste. Además, la flexibilidad del algoritmo se mejora al dar la oportunidad al usuario de incluir algunos parámetros autodefinidos, como la función de pérdida, o la métrica utilizada para la validación y prueba (^{Urraca, Antonanzas, Antonanzas-Torres, & Martinez-De-Pison, 2017}).

Implementación del modelo XGBoost en la estimación de la ETo

Para la implementación del modelo XGBoost, como primer paso se optimizaron los hiperparámetros nrounds, max_depth, eta, gamma, colsample bytree, min_child_weight y subsample (Tabla 4) utilizando la librería Caret del Software R (^{Fan et al., 2018a}); segundo, se ajustó el modelo XGBoost utilizando la librería “Xgboost” del Software R, utilizando validación cruzada (VC = 5 fólders) para evitar el sobreajuste. El 70 % de los datos se usó para la etapa de entrenamiento y el 30 % para la validación. Las variables de entrada al modelo son los valores de los datos meteorológicos de T_máx, T_mín y R_a , los valores de ETo-FAO56PM como variable objetivo.

Análisis estadístico

En el presente estudio se utilizaron cuatro indicadores estadísticos para evaluar el desempeño de los modelos convencionales y de inteligencia artificial, estos indicadores son coeficiente de determinación (R²; Ecuación (13)); raíz cuadrada media del error (RMSE; Ecuación (14)); error absoluto medio (MAE; Ecuación (15)); error medio de sesgo (MBE; Ecuación (16)):

donde n es el número de comparaciones; P_i y O_i , valores estimados y observados de ETo-FAO56PM, respectivamente; P_prom , el promedio de los valores estimados de ETo; O_prom , el promedio de los valores observados de ETo; las unidades de ETo se encuentran en mm d^-1.

R² se usa por lo común para estimar el desempeño de modelos hidrológicos. Representa la fracción de los valores estimados que son los más cercanos a la línea de datos de medición. Valores del coeficiente de determinación cercanos a 1 indican modelos más eficientes y la línea de regresión se ajusta mejor a los datos. EL RMSE es una medida utilizada con frecuencia para comparar errores de predicción en diferentes modelos; cuanto menor sea su valor, mejor será la capacidad predictiva de un modelo en términos de su desviación absoluta. El MAE es la suma de los valores absolutos de los errores divididos por el número de observaciones; se usa con frecuencia para medir qué tan cerca están los valores estimados de los valores observados. El MBE proporciona información sobre la tendencia del modelo a sobreestimar o subestimar la variable; cuantifica el error sistemático del modelo.

La evaluación de la incertidumbre de los indicadores estadísticos R², RMSE, MAE y MBE se realizó mediante la construcción de intervalos de confianza bootstrap (ICB) al 95 % de nivel de confianza. Para tal propósito se utilizó el método no paramétrico bootstrap percentil como técnica de remuestreo, utilizando B = 10 000 réplicas con remplazamiento, con el propósito de inducir mayor precisión en las estimaciones (^{Efron, 1992}).

Los ICB ofrecen una manera de estimar con alta probabilidad un rango de valores en el que se encuentra el valor del parámetro (indicador estadístico).

También se computó el error estándar de la distribución (ee) y el valor de cada indicador estadístico bootstrap, calculando la desviación estándar y la media de las B réplicas.