Introducción

La precipitación es un fenómeno caótico que fluctúa en tiempo y espacio, por lo que conocer el comportamiento de sus patrones es una tarea sumamente complicada, pero necesaria para los administradores de los recursos hídricos debido a que existe la necesidad de prevenir o planificar eventos normales o extremos como sequías o inundaciones. Las series de tiempo de registros hidroclimatológicos tienen la característica esencial de ser no estacionarias y oscilatorias, o no lineales. Sin embargo, algunas veces por necesidad y para dar solución a un problema específico, se tratan como estacionarios y con tendencias lineales. Los modelos usados de manera tradicional han sido descritos en la literatura, como regresión lineal múltiple (RLM) y autorregresivos integrados de media móvil (ARIMA) (^{Adamowski, 2008}; ^{Dastorani, Afkhami, Sharifidarani, & Dastorani, 2010}; ^{Nourani, Alami, & Aminfar, 2009}; ^{Remesan, Shamim, Han, & Mathew, 2009}).

En las últimas décadas se han importado de otras especialidades nuevas tecnologías de análisis de señales y pronóstico, como las Artificial Neural Networks (ANN) y sus variantes. En estas redes no se involucran las variables físicas implicadas en los fenómenos, pero se pueden representar la no linealidad y no estacionalidad de los datos meteorológicos que desempeñan un papel significativo en el pronóstico (^{Adamowski, 2008}). Las ANN matemáticamente son a menudo vistas como aproximadores universales (^{Govindaraju, 2000}) y son caracterizadas por su arquitectura (que representa el patrón de conexión entre los nodos), por su método para determinar el enlace entre los pesos, y la función de activación (^{Govindaraju, 2000}). Las ANN son muy hábiles para encontrar un patrón y predecirlo, además han sido utilizadas de manera extensa para la predicción de diversas variables en el área de recursos del agua: caudales, precipitación, sedimentos suspendidos, modelos lluvia-escurrimiento, etcétera (^{Besaw, Rizzo, Bierman, & Hackett, 2010}; ^{Dastorani et al., 2010}; ^{Kişi, 2009}; ^{Partal, 2009}).

^{Jang (1993)} introduce Adaptive-Networkbased Fuzzy Inference System (ANFIS). Este método es una adaptación de ANN con Fuzzy Interfase System (FIS) para el modelado de los aspectos del conocimiento humano y procesos de razonamiento sin utilizar análisis cualitativo preciso, con el propósito de que se puedan interpretar de manera sencilla con reglas difusas del tipo if-then. Un FIS emplea funciones de membresía (MF), por ejemplo μA(x), que corresponde al grado en el cual el valor de x tiene el atributo A (pequeño, largo, etc.). Los parámetros del FIS están determinados por los algoritmos de aprendizaje de la ANN. La parte premisa es la entrada del modelo y la consecuente es la salida del modelo. Existen tres tipos de ANFIS. En el tipo I, la salida global es el promedio ponderado de las salidas de cada regla inducida por la fuerza de disparo de la regla y la salida de la MF; en el tipo II, la salida difusa global se deriva por la aplicación de una operación "máx" a la salida difusa calificada; en el tipo III se usan reglas difusas if-then Takagui y Sugenno, la salida de cada regla es una combinación lineal de las variables de entrada más una constante y la salida final es el promedio ponderado de la salida de cada regla.

Recientemente, ^{Besaw et al. (2010)} predijeron el caudal de una cuenca no aforada. ^{Wu, Chau y Fan (2010)} usaron ANN para predecir series de tiempo de lluvia con datos mensuales y diarios. ^{Dastorani et al. (2010)} aplicaron varias técnicas de ANN y ANFIS para la predicción de la precipitación registrada por una estación en el centro de Irán con clima hiperárido. ^{Silverman y Dracup (2000)} también emplearon ANN para predecir los patrones de lluvia en California, EUA.

En México se han utilizado las ANN en problemas hidrológicos. ^{Molina-Aguilar y Aparicio (2006)} utilizaron las FFBP para simular el tránsito de avenidas. Por su parte, ^{Cisneros-Iturbe y Pelczer (2007)} utilizaron las ANN para pronosticar alturas de precipitación. ^{González-Camacho, Cervantes-Osornio, OjedaBustamente y López-Cruz (2008)} predijeron la evapotranspiración de referencia de los registros históricos elaborados en el estado de Sinaloa.

Sin embargo, es importante tener en cuenta lo mencionado en ASCE (2000b): la precipitación sirve como motor para la mayoría de los procesos hidrológicos, además predecirla es una tarea complicada, porque presenta un grado amplio de variabilidad espacial y temporal. ^{Partal y Kişi (2007)}, ^{Kişi (2009)}, y ^{Cisneros-Iturbe y Pelczer (2007)} mencionan que la predicción exacta de la precipitación es muy difícil debido a la complejidad de los procesos físicos involucrados. Por un lado, Partal y Kişi (2007) indican que lo anterior ocurre en especial al comparar los datos de escorrentía; por su parte, ^{Cisneros-Iturbe y Pelczer (2007)} refieren que la correlación entre los valores de precipitación sucesivas puede ser menor en relación con el mismo análisis de los valores de escorrentía; por ello, el número de publicaciones para el pronóstico de precipitaciones es escaso.

La transformada wavelet (TW) fue introducida por ^{Grossmann y Morlet (1984)} como una herramienta de análisis tiempo-escala para señales no estacionarias, así como para estudiar las no linealidades y discontinuidades de datos irregulares. También puede usarse para analizar series de tiempo que contienen no estacionalidades en muchas frecuencias diferentes (^{Daubechies, 1990}). La TW consiste en descomponer una serie de tiempo en el espacio tiempo-frecuencia y determinar sus modos dominantes de variabilidad, y cómo éstos se modifican en el tiempo (^{Torrence & Compo, 1998}). En este proceso de transformación se mantiene la ubicación de la señal en el tiempo y es posible reconstruirla (^{Graps, 1995}; ^{Percival & Walden, 2006}). La TW compara una señal con versiones dilatadas y trasladadas de ciertas funciones wavelet llamadas wavelet base o madre (^{Graps 1995}; ^{Remesan et al., 2009}) que deben satisfacer ciertos requisitos matemáticos (^{Graps, 1995}). Existe una variedad de estas funciones, entre ellas Haar, Morlet, Daubechies, Symlets, Meyer y Coiflets (ver figura 1).

a) Haar, b) Morlet, c) Daubechies.

Figura 1 Diferentes tipos de funciones wavelet base o madre. 

La teoría wavelet se divide en dos tipos de transformaciones: TW continua (TWC) y TW discreta (TWD) (^{Addison, Murray, & Watson, 2001}; ^{Daubechies, 1992}). La TWC es en esencia una herramienta de análisis exploratorio de datos que puede ayudar al ojo humano a distinguir características de su interés (^{Percival & Walden, 2006}). Para explicar el efecto de borde, este trabajo utiliza la TWD.

De forma reciente han surgido nuevas técnicas de preprocesamiento para ayudar en el pronóstico de algunos eventos climáticos, como la precipitación y sus variables asociadas, atendiendo el conocimiento profundo de la no linealidad del sistema climático, por ejemplo ^{Besaw et al. (2010)}, ^{Remesan et al. (2009)} y ^{Wu et al. (2010)}. Esta nueva perspectiva es principalmente desarrollada con la TWD, en conjunto con una herramienta estadística; ambas son utilizadas como métodos de preprocesamiento de datos. Este preprocesado, acoplado a algún algoritmo de pronóstico, principalmente ANN o cualquiera de sus variantes, ha ganado terreno en el pronóstico de series climáticas sobre otros métodos (^{Adamowski & Sun, 2010}; ^{Kişi, 2009}; ^{Partal, 2009}; Remesan et al., 2009; ^{Shiri & Kisi, 2010}). El acoplamiento de ambos procesos se conoce como modelo wavelet-neuronal (WN).

Existen tres enfoques de preprocesamiento reportados en la literatura. El primero se basa en descomponer la señal con TWD, y utilizar las descomposiciones como insumos a las ANN o alguna de sus variaciones (^{Kişi, 2009}; ^{Nourani et al., 2009}; ^{Partal, 2009}; ^{Remesan et al., 2009}). El segundo enfoque consiste en suavizar la señal y una vez encontrada una nueva señal se predice su comportamiento, siendo ejemplos los mencionados por ^{Partal y Kişi (2007)}, ^{Shiri y Kisi (2010)}, ^{Partal (2009)}, ^{Wu, Chau y Li (2009)}, y ^{Tiwari y Chatterjee (2010)}. Este enfoque tiene el problema de pérdida de información en el proceso del suavizado. ^{Remesan et al. (2009)} concluyen que es más conveniente trabajar con los datos crudos debido a que se trabaja siempre con la señal original sin modificarla o suavizarla. El tercero es de descomposición-pronóstico-reconstrucción. Se descomponen las series en subseñales con la TWD. Cada una de ellas se predice utilizando alguna herramienta para este fin y, por último, a partir de las subseries pronosticadas, se reconstruye la serie original, utilizando la transformación inversa o ANN; esta última técnica la utilizan ^{Kim y Valdés (2003)}. Los enfoques presentan el problema del fenómeno del efecto de borde producto de la transformada. El inconveniente es que precisamente estos datos del final de la serie son los que se utilizan para realizar el pronóstico. Pese a la importancia del fenómeno para la predicción, no se menciona en la literatura. Por ejemplo, ^{Kim y Valdés (2003)} dicen que utilizan un método de espejo, y ^{Adamowski (2008)} indica un método para el corte de datos y bordes. A fin de explicar este fenómeno es necesario citar el método de la convolución, procedimiento que utiliza la transformada para llevarse a cabo. En ella, los coeficientes wavelet utilizados para el filtrado son trasladados a lo largo de la señal. Al inicio se emplean valores que se sitúan adelante de la serie, esto es, los futuros; al final, la señal se acaba, pero el proceso aún no ha terminado; para completarlo se requiere obtener una serie de valores; para solucionarlo se rellenan de forma arbitraria esos espacios vacíos mediante algún procedimiento sencillo.

En este trabajo se analizan los dos primeros enfoques de pronóstico y se desarrolla con los algoritmos de la MODWT y de Mallat, aplicando FFBP y ANFIS a datos de precipitación mensuales, a fin de realizar un pronóstico inmediato de la precipitación con un mes de anticipación en una región con clima seco templado del estado de Baja California, México. Se discute la influencia del efecto de borde en el pronóstico.

Metodología

Se ha escogido la información de la Estación Climatológica 016 Ensenada (latitud 31° 53' 45" y longitud 116° 35' 47"), que se encuentra al noroeste del estado de Baja California Norte, en la vertiente del océano Pacífico, situada en la cortina de la presa Emilio López Zamora de la ciudad de Ensenada, que en lo sucesivo se nombrará como Ensenada (figura 2). De ella se escogieron los registros de precipitación recopilados por la Comisión Nacional del Agua (Conagua). La estación se encuentra en una región de clima seco templado, con un periodo de lluvias de octubre a marzo (^{INEGI, 1995}). La serie de tiempo de precipitación se elaboró con los registros de mayo de 1948 a febrero de 2011, a intervalos de medición mensual, lo que corresponde a 754 valores de los cuales se utilizan 622 para entrenamiento, los siguientes 100 para validación y los restantes 32 para pronóstico.

Figura 2 Localización de las estaciones Abelardo L. Rodríguez en Tijuana y Emilio López Zamora en Ensenada, Baja California Norte. 

FFBP

El FFBP, ampliamente utilizado en hidrología (^{Adamowski, 2008}; ^{Cisneros-Iturbe & Pelczer, 2007}; ^{González-Camacho., 2008}; ^{Kim & Valdés, 2003}; ^{Kişi, 2009}; ^{Nourani., 2009}; ^{Partal, 2009}), es una entrada y dos capas. Proporciona un mapeo funcional no lineal entre el grupo de variables de entrada y salida. Está basado en una combinación lineal de variables de entrada, que son transformadas por una función de activación no lineal (^{Kim & Valdés, 2003}).

La expresión explícita de un valor de salida de una FFBP es dada por:

(1)

donde wji = peso de conexión de la capa oculta de la i-ésima neurona de la entrada a la j-ésima neurona de la capa oculta; wjo = sesgo de la j-ésima neurona oculta; fh = función de activación de la capa oculta; wkj = peso de conexión de la capa de salida de la j-ésima neurona de la capa oculta a la k-ésima neurona de la capa de salida; wko = sesgo para la k-ésima neurona de salida, y fo = función de activación de la capa de salida. Los pesos son diferentes para cada conexión y sus valores pueden ser cambiados durante el proceso de entrenamiento (^{Kim & Valdés, 2003}). En la figura 3 se presenta de modo esquemático una ANN FFBP alimentada hacia adelante con entrenamiento de retropropagación de los errores.

Figura 3 ANN FFBP alimentada hacia adelante con entrenamiento de retropropagación de los errores. 

Se utilizó una entrada, una capa oculta y una de salida. Una sola capa oculta es suficiente para alcanzar un buen desempeño de la FFBP (^{Beale & Jackson, 1990}; ^{Cisneros-Iturbe & Pelczer, 2007}), y es determinada por el número de neuronas y sus conexiones (^{Beale & Jackson, 1990}). La cantidad de estas neuronas de la capa oculta fue definida por ensayo y comprobación, como se hace por lo general (^{Govindaraju, 2000}; ^{Kişi, 2009}). Asimismo, se empleó el algoritmo de Levenberg-Marquardt como función de entrenamiento debido a que tal técnica es más poderosa y rápida en este tipo de funciones (^{Braham & Varas, 2003}; ^{Kişi, 2009}). Además, ha sido utilizada por ^{Partal y Kişi (2007)}, ^{Remesan. (2009)}, ^{Tiwari y Chatterjee (2010)}, y ^{Wu. (2010)}, entre otros.

En el presente estudio, el número de neuronas de entrada es variable y depende del mes actual y los meses de antelación considerados para cada configuración; mientras que el número de neuronas de la salida es siempre uno, es decir, el mes posterior al mes actual. Las entradas y salidas se normalizan en el rango de0.1 a 0.9 antes de introducirlo en la red.

ANFIS

Para explicar su funcionamiento se debe consultar la figura 4. Por simplicidad se asume un FIS con dos entradas x y y, y una salida f. Se supone que la regla base contiene dos reglas difusas if- then del tipo Takagi y Sugeno. Las funciones de los nodos en la misma capa son iguales:

Figura 4 a) Modelo difuso de primer orden tipo Sugeno, con dos estradas y dos reglas; b) arquitectura equivalente de ANFIS. 

(2)

donde x y y son las entradas al nodo i (i = 1, 2), Ai y Bi son niveles lingüísticos caracterizados por una función de membresía apropiada; pi, qi y ri son los coeficientes de la combinación lineal en el sistema de inferencia tipo Sugeno; también son llamados parámetros consecuentes.

Capa 1. Nodos de entrada, para cada nodo existe una función de membresía (ecuación (3)), que arroja un valor lingüístico:

(3)

donde0i1 es la MF de Ai y especifica el grado en que x satisface el cuantificador Ai, el superíndice uno significa la capa. Los parámetros ai, bi y ci (parámetros premisa) son los parámetros de la función de la MF, que cambian su forma; por lo usual se utiliza una campana de Gauss con valores de 1 y 0 como máximo y mínimo, respectivamente.

Capa 2. Reglas de los nodos. Un operador and/or es aplicado para obtener una salida que es llamada (0i2) fuerza de disparo de una regla (ecuación 4):

(4)

Capa 3. Promedio de los nodos. Se calcula la relación de fuerza de disparo de la i-ésima regla con respecto a la suma de las fuerzas de disparo de todas las reglas y las salidas de esta capa son llamadas fuerzas de disparo normalizadas 0i3:

(5)

Capa 4. Nodos consecuentes:

(6)

0i4 es la contribución de la i-ésima regla en relación con el total de salidas del modelo.

Capa 5. Nodo de salida. Es una capa de un solo nodo y se calcula la salida final como la suma de todas las señales entrantes 0i5:

(7)

ANFIS utiliza un algoritmo de aprendizaje híbrido. Cuando se encuentra en el recorrido hacia adelante, las señales siguen hasta la capa 4 y los parámetros consecuentes son estimados por el método de los mínimos cuadrados. En el recorrido hacia atrás, las tasas de error se propagan en ese sentido y los parámetros premisa son actualizados por el método del gradiente descendente. Se puede encontrar más información en ^{Jang (1993)}.

Transformada Wavelet

La TW es una descomposición tiempo frecuencia, que usa funciones wavelet localizadas en tiempo y frecuencia (^{Daubechies, 1990}). La TWC de una función s(t) es definida como:

(8)

donde ψ(t - b)/a es una función wavelet admisible y se utiliza para analizar la señal s(t). Los coeficientes de transformación T(a, b) son encontrados para una específica localización, t = b, y para un específico periodo wavelet, que es función de a, el parámetro de escala. El diagrama de T(a, b), con sus respectivos parámetros a y b es conocido como escalograma. Como en la transformada de Fourier, la señal original puede ser reconstruida usando la inversa:

(9)

donde Cg es conocida como constante de admisibilidad.

^{Mallat (1999)} propuso discretizar los valores de a y b, asignando 2 y 1, respectivamente. Con esta modificación construyó la base ortogonal wavelet. La diádica se escribe como sigue:

(10)

donde j y n son enteros que modifican la función wavelet; el primero es el parámetro de escala que controla la dilatación y el segundo la translación. T_m,n es el coeficiente wavelet en esa escala y posición. Si si es una serie de tiempo de longitud N potencia de 2, el subíndice i = 0, 1, 2,... N - 1; si ocurre en un tiempo discreto i. La TWD inversa se puede obtener mediante la siguiente ecuación:

(11)

donde T es la señal media en el caso de N = 2^j

El concepto multirresolución consiste en aplicar sistemáticamente la TW. Para el caso práctico se utilizan filtros: la función de escalamiento φ(x) y la función wavelet ψ(x) que son relacionadas, respectivamente, a los filtros pasa-bajas h y pasa-altas g, que en realidad los dos son pasa-bandas (^{Percival & Walden, 2006}):

(12)

(13)

donde t = 0, 1, ..., NJ - 1, A_m,t es la aproximación de la señal en la escala m, posición diádica t, D_m,t es el detalle para los mismos parámetros. L es la longitud del filtro y el índice (m - 1, 2t + 1 - l mod N_m-1 ) significa el submuestreo y remuneración diádica, y la función mod es el módulo aritmético para realizar la convolución. El proceso de reconstrucción es la transformada inversa.

La MODWT no es una transformada ortonormal. El procedimiento para llevarla a cabo es similar al desarrollado por el algoritmo de Mallat. Se puede encontrar más información en ^{Percival y Walden (2006)}. Las descomposiciones se realizan de la siguiente manera:

(14)

(15)

donde t = 0, 1, ..., N -1. El símbolo à es la aproximación del algoritmo de la MODWT; D es el detalle;

La señal original puede ser recuperada usando la transformación inversa.

Convolución

La convolución es una operación matemática formal, tal como la multiplicación, suma e integración, y es utilizada en muchos campos de las matemáticas; además, se emplea en el procesamiento digital de señales debido a que relaciona la señal de entrada, la señal de salida, y la respuesta o impulso. La convolución discreta es el proceso de invertir una secuencia, desplazarla un paso a la derecha y sumar los productos superposicionados; también es llamada multiplicación en serie. Para explicar este documento es importante regresar a las raíces de este proceso debido a que la DWT tiene implícita una operación de este tipo.

En la figura 5 se muestra de modo esquemático el desarrollo de la convolución para la obtención de las aproximaciones con el algoritmo de DWT Mallat; esto es, se utiliza el filtro g con cuatro coeficientes de escalamiento c_k . La secuencia de los cuatro dígitos 1/√2c_k, filtra la señal S_0,n en escala m = 0 en submuestreos en cada salto; se encuentra una componente de aproximación en la siguiente escala m + 1 y se desliza el filtro pasa bajas utilizando dos saltos en la posición de la aproximación para encontrar la aproximación A_1,n' en escala m = 1. El submuestreo remueve todos los segundos valores; el esquema muestra los valores retenidos en la escala m = 1. Para los fines de este trabajo es importante observar que en cualquier posición t se utilizan tres valores siguientes en cada salto para encontrar un valor en la siguiente escala m + 1; es decir, en el proceso de convolución, los coeficientes de escalamiento y wavelet utilizados para el filtrado son trasladados a lo largo de la señal; en cualquier posición t se emplean los valores que se sitúan adelante de la serie o sea los futuros de esa posición. En la figura 5, la señal de ocho elementos (S _0.0, S _0.1,..., S _0.7) en el nivel m = 0 es filtrada por cuatro coeficientes de escalamiento (c ₀, c ₁, c ₂, c ₃). Al final, para terminar este proceso, es decir, para encontrar el elemento A _1,3' la señal se acaba, pero el proceso aún no ha terminado; entonces, para completarlo, se requiere obtener una serie de valores, es decir, se rellena la serie. En nuestro ejemplo, los coeficientes de escalamiento se multiplican por S _0.6' S _0.7' l _0.1' l _0.2' los últimos dos elementos l _0,X que aparecen se colocan de manera arbitraria; en otras palabras, se infieren de acuerdo con el criterio del investigador. Para ello existen varios procedimientos clásicos. El primero consiste en realizar un espejo de los últimos valores, como lo utilizan ^{Kim y Valdés (2003)}, y ^{Adamowski (2008)}. El segundo se basa en repetir el inicio de la serie. El tercero es asignar ceros. El cuarto se coloca en la media de la señal. En fin, se pueden poner los valores que se deseen. Esto tiene una implicación en el pronóstico debido a que se convierte en una predicción implícita sobre la que se desea realizar, pero en este caso se colocan valores arbitrarios y en el otro se emplea un método robusto de predicción como la FFBP o ANFIS. En relación con otros algoritmos para realizar la transformada, el "a trous" presenta mayor penetración en la señal, específicamente es el número de coeficientes por 2 2 por cada nivel; en el MODWT, en cada nivel se introduce el número de coeficientes menos uno; y el de Mallat, conforme se baja de nivel avanza el número de coeficientes de escalamiento menos dos; pero producto de las decimaciones, la señal se va reduciendo en mitades cada vez. Este último algoritmo con el wavelet de Haar no presenta efecto de borde debido a que sólo tiene dos coeficientes de escalamiento y wavelet (^{Cornish, Bretherton & Percival, 2006}). En el proceso inverso de reconstrucción del algoritmo de Mallat, el filtrado es simplemente invertido, por lo cual los componentes de la escala más grande son regresados a través de los filtros. Para el caso de la MODWT, el método es similar al de Mallat, la diferencia es simplemente que se omite el submuestreo en el proceso de filtrado.

Figura 5 Proceso de convolución en el filtrado de descomposición de la aproximación con el algoritmo de Mallat. 

Modelo wavelet-neuronal/wavelet neurodifuso

Primer enfoque

Una vez obtenida la serie de tiempo de la precipitación mensual de la estación en estudio, se divide en tres grupos de valores: entrenamiento, validación y pronóstico. El primero se utiliza para entrenar el algoritmo de pronóstico; el segundo se emplea para validar o verificar en el proceso de entrenamiento de un grupo de datos para los cuales la red no se está calibrando, y el tercero son los valores que se utilizan para verificar el pronóstico. A continuación se descomponen las series hasta el nivel m del 1 al 7 mediante el clásico algoritmo de Mallat o MODWT con las funciones wavelet de Haar, Daubechies de orden cuatro (Db4); estas descomposiciones se utilizan como insumos al FFBP o ANFIS, según sea el caso; en específico, se introducen los valores de antelación de la señal, la última aproximación de ésta y los detalles. Para demostrar la influencia del efecto de borde se utilizaron cinco técnicas. La primera consistió en descomponer toda la serie, con el propósito de que este fenómeno no influyera en el resultado y se denomina sin efecto de borde (SEB); pese a ello, esto no es una predicción debido a que se está incluyendo la información del periodo de predicción y siempre es necesario tener sumo cuidado en no incluir información futura (^{Kim & Valdés, 2003}). La repetición de los datos del inicio de la serie fue la segunda y se basa en colocar después del último valor de la serie, el primero, segundo y así de manera sucesiva, hasta completar el número de datos requeridos por el proceso. En tercer lugar se utilizó un espejo de los valores, esto es, poner después del último valor de la serie el último, penúltimo, antepenúltimo, y si sucesivamente, hasta lograr la transformada. En la cuarta se colocan los valores de la media en esos valores. Por último, la quinta técnica consistió en ubicar al final de la serie los valores 422, 423, 424, etcétera, debido a que ese grupo de valores tiene el coeficiente de correlación más alta (0.823) con el grupo de datos de pronóstico. Un esquema de tal enfoque se encuentra en la figura 6.

Se descompone la señal S(t) mediante la DWT, los componentes de detalle SD1(t), SD2(t),..., SDj(t); la última aproximación SAj(t); y la señal original S(t), con diferentes valores de antelación, son utilizados como entradas a la ANN o ANFIS, y mediante los pesos W se encuentra el siguiente valor de la serie S(t + 1) pronosticado.

Figura 6 Esquema del primer enfoque. 

El proceso de pronóstico consiste en descomponer la serie original hasta el valor t. Se aplica el algoritmo de pronóstico y se obtiene el valor siguiente t + 1 pronosticado; la segunda vez se descompone hasta el valor t + 1 de la serie original y se obtiene el siguiente valor t + 2 pronosticado, y así de forma sucesiva hasta llegar al final de la serie original. Es importante recalcar que el pronóstico es a plazo inmediato, o sea, se encuentra el valor siguiente. De acuerdo con la literatura revisada, la precipitación es complicada de predecir con estas técnicas a corto, mediano o largo plazos.

Para encontrar las configuraciones óptimas se realizó un mapeo multifuncional. Para el caso del FFBP, se llevaron a cabo pruebas con diferentes configuraciones (1, 0.5, 1.0, 1.5 n) (^{Kim & Valdés, 2003}) para definir el número de neuronas de la capa oculta, donde n es igual al número de neuronas de la entrada. Además, se utilizó el algoritmo de terminación temprana (^{Cisneros-Iturbe & Pelczer, 2007}; ^{Govindaraju, 2000}) para hacer un pronóstico óptimo. El máximo número de épocas fue 250. Los valores de entrada a la red fueron normalizados de 0.1 a 0.9 para ser usados por la red y se repitió el entrenamiento de la red hasta que arrojara un resultado consistente debido a que los pesos se inicializan de modo aleatorio. Para el caso del ANFIS, se varió el número de funciones de membresía de 2 a 6, que consistieron en una campana de Gauss; además, se varió el número de meses de antelación de 1 a 6, siendo uno el mes actual; se utilizó un vector de validación y se modificó el número de épocas de 1 a 25. Es común utilizar como medida de desempeño de un algoritmo de pronóstico la Root Mean Square Error (RMSE), en este trabajo se emplea y su valor objetivo varió de 1x10^-1 a 1x10^-9, con intervalos de la unidad en escala logarítmica; además, se utilizó el coeficiente de determinación R ² para ambos casos.

Segundo enfoque

Este método consiste en modificar la señal, con el propósito de ayudar al algoritmo de predicción a modelar una serie más suave; para ello se eliminan las altas frecuencias. Para lograrlo se descompone la señal normalmente en etapa de reconstrucción. Se eliminan de forma sistemática los detalles del nivel uno hasta el m elegido, lo que obtiene una nueva señal más suave. Al eliminar esos detalles se pierde una parte de la señal; se quita una fracción del efecto de borde y produce uno nuevo; además, esa fracción que se excluyó produce un ruido que tiene la forma del wavelet madre; este fenómeno es más significativo en relación con el número de detalles que se eliminen y su magnitud. Una vez que se obtiene la nueva señal, los algoritmos de predicción trabajan como si fueran autónomos. Para este caso, se utilizan los mismos criterios que el enfoque anterior, con las mismas configuraciones; de igual manera, se trabaja para la TWD, pero en esta ocasión se utiliza sólo el algoritmo de Mallat con la función wavelet Db4. Los criterios de desempeño se obtienen de acuerdo con la señal original, que es el objeto de nuestra predicción. Con el propósito es identificar cuáles coeficientes wavelet de descomposición presentan mayor similitud con la señal original, se calculó el coeficiente de correlación para las componentes de los niveles del 1 al 9. La figura 7 muestra un esquema de este enfoque.

Se descompone la señal S(t) mediante la DWT; después se eliminan de forma sistemática los coeficientes de detalle; se reconstruye la señal mediante la DWT inversa, de esta manera se encuentra una señal nueva más suave S _Nueva(t), que se utiliza como entrada con diferentes valores de antelación a los algoritmos de ANN o ANFIS, y mediante los pesos WS _Nuevos se encuentra el siguiente valor de la serie nueva S_N (t + 1) pronosticado.

Figura 7 Esquema del segundo enfoque. 

Resultados y discusión

En el cuadro 1 se presentan los valores estadísticos de las series de tiempo de la estación Ensenada. La primera columna muestra la estación con cuatro valores: número de datos de entrenamiento, validación, pronóstico y el total de ellos. De la 2 a la 11 se exhiben para cada grupo de datos la media, valor máximo, mínimo, desviación estándar, sesgo, coeficiente de autocorrelación 1, 2, 3, 12 y 24 meses, respectivamente, y por último el exponente de Hurst.

Cuadro 1 Datos de precipitación de las estaciones Ensenada y Tijuana. 

En el cuadro 1, la media del total de los datos es de 21.71 mm; sin embargo, la media para los grupos de datos de entrenamiento, validación y pronóstico son 21.47, 18.64 y 27.09, respectivamente. Esto significa que en el proceso de entrenamiento de la FFBP o ANFIS se utilizará un grupo de valores que tienen una media (21.47) mayor que los datos de validación (18.64); a pesar de ello, en la etapa de predicción, el grupo de valores de pronóstico tiene una media de 27.71 mayor que los datos de entrenamiento. Por otra parte, el valor máximo de precipitación se encuentra en el grupo de datos de entrenamiento con 237.6 mm, el valor intermedio para los datos de validación con 171.3 mm y el más pequeño en los de pronóstico con 145.5 mm Todas las distribuciones están sesgadas hacia la derecha, aunque las más sesgadas son las de entrenamiento. La autocorrelación uno es la más significativa en relación con las más cercanas al final de la serie, y las autocorrelaciones 12 y 14 indican una cierta periodicidad anual; es decir, existe similitud con los valores anuales. Estos resultados preliminares muestran que los algoritmos de pronóstico trabajarán los tres grupos de valores, que son heterogéneos desde el punto de vista de la estadística y el exponente de Hurst, que es el único indicador no lineal o caótica indica que la serie de tiempo Browniana o al azar; es decir, no se puede predecir debido a que su valor es cercano a 0.5 en todos los casos.

En el cuadro 2 se presentan los resultados del modelo del primer enfoque. La primera columna es el modelo de predicción con el método de rellenado en el borde. La segunda es la función wavelet. La tercera se refiere al nivel de descomposición de la transformada wavelet. La cuarta es el resultado en términos del RMSE del mejor entrenamiento. La quinta es la arquitectura del modelo de predicción. La sexta es el error objetivo. La séptima son las épocas o iteraciones empleadas. De la octava a la última se muestran los indicadores de desempeño RMSE y R ² para los tres grupos de valores entrenamiento, validación y pronóstico. Se presentan los mejores resultados de los mapeos realizados con todas las configuraciones y los más simples debido a que existen resultados similares con otros arreglos.

Cuadro 2 Resultados del modelo WN del primer enfoque. 

En la primer columna se muestran los siguientes algoritmos FFBP; a continuación, ANFIS, ambos autónomos; después, la media aritmética. Además, se muestran tres grupos de valores Mallat-FFBP, que significa la descomposición mediante la TW con el algoritmo de Mallat y la predicción con FFBP. El siguiente grupo, Mallat-ANFIS, es similar al anterior, pero el algoritmo de predicción con ANFIS y estos dos grupos anteriores utilizan cuatro técnicas de rellenado SEB. En este caso, se rellena la serie con los valores futuros, es decir, el efecto de borde no influye. Repetir, espejo, media son como se explicó en la convolución. Los últimos dos grupos difieren de los anteriores en la manera de realizar la transformada, pues estos últimos se descomponen empleando la MODWT. Para demostrar la influencia del efecto de borde se añade una técnica denominada 422 debido a que desde este valor el grupo de pronóstico tiene un mayor coeficiente de correlación.

La función wavelet Db4 para los casos de la configuración SEB es la mejor opción porque que la penetración del efecto de borde guarda la información y le ayuda al algoritmo de pronóstico a conocer la evolución futura. La función wavelet de Haar no presenta efecto de borde, pues con la descomposición de Mallat y al tener dos coeficientes de escalamiento y wavelet no requiere de rellenado. Es la mejor opción para la configuración Mallat-FFBP con las opciones repetir, espejo y media; de hecho, si se compara con el ANFIS autónomo hay resultados similares para los grupos de valores de entrenamiento y pronóstico. La validación entrega resultados erráticos. Para la mayoría de los casos, las descomposiciones wavelet son para niveles 1 y 2, excepto para la configuración MODWT-FFBP, que son en los niveles del 3 al 5.

En todos los casos existe un mejor entrenamiento en términos del RMSE, pero éste muestra resultados en el pronóstico más pobres que los presentados en la columna ocho, por esta razón se ajustó el error objetivo, para que la red no quedara apretada y perdiera generalización en el proceso de pronóstico. Por otra parte, al utilizar en la FFBP la función de entrenamiento de Levenberg-Marquardt, se redujo de modo significativo el número de épocas. Para el caso de ANFIS, la cantidad de épocas fue de diez para todos los casos.

En ese mismo cuadro 2 nos referiremos ahora a los métodos híbridos, en específico a las configuraciones repetir, espejo y media; en todos los casos, los resultados de pronóstico en términos de RMSE son superiores a 30 mm, pero en la etapa de entrenamiento algoritmo de Mallat, entrega resultados inferiores que los de MODWT. Esto se debe a que este último algoritmo tiene mayor penetración del efecto de borde, entonces existe mayor información futura almacenada en los coeficientes que se introdujo en el momento de realizar la convolución. En la configuración MODWT-FFBP con el rellenado de espejo entrenamiento perfecto en términos de la pareja (RMSE, R2) (1.11, 1.0) , validación (0.12,1.0), pero en pronóstico (40.93, 0.08), este grupo de datos tiene un coeficiente de correlación de 0.19; sin embargo, se utiliza otro grupo de datos, por ejemplo, el 422 en el rellenado con un coeficiente de correlación mayor (0.823); la pareja de indicadores de desempeño se incrementa a (23.68, 0.67). Ahora, si utilizamos en el rellenado un grupo de valores con coeficiente de correlación igual a 1.00, es decir, la configuración SEB, la pareja de indicadores de desempeño es de 0.00, 1.00. Dicho de otra forma, se puede usar un rellenado de espejo con una serie más periódica de tal forma que el coeficiente de correlación se incremente y automáticamente el indicador de desempeño también suba; es decir, los métodos híbridos y el pronóstico mostrados están íntimamente relacionados con el coeficiente de correlación que existe entre los valores de relleno y los de pronóstico debido a que esta información está guardada en los coeficientes del proceso de convolución. De todos los métodos de pronóstico válidos, esto es, que no introducen información futura, los resultados son coherentes con el indicador caótico de Hurst; es decir, la serie de tiempo no se puede predecir, es al azar; por lo tanto, los indicadores de desempeño exponen resultados bajos, como era de esperarse.

En la figura 8 se presentan los resultados del primer enfoque para dos casos empleando FFBP. El primero es con el rellenado (SEB) y el pronóstico es prácticamente perfecto (6.39, 0.97), evidente; el coeficiente de correlación es 1.00, y para el segundo caso se utiliza MODWT, con el rellenado espejo; en este caso, se emula el comportamiento, pero con cierto desfase (40.93, 0.08). Esto se debe a que los datos de relleno están desfasados, es decir, el coeficiente de correlación es 0.18. En la figura 9, todos los casos son con ANFIS, el primero es cuando es autónomo y el segundo cuando se utiliza MODWT con el rellenado (422) y el pronóstico presenta una pareja de valores (23.68, 0.67); esto es, se incrementa de forma significativa el índice de desempeño. Lo que sucede es que se elimina el desfase debido a que el grupo de valores de relleno están más correlacionados (0.826) con los de pronóstico.

Figura 8 Resultados en la etapa de pronóstico los modelos Mallat-FFBP (SEB) y MODWT-FFBP (espejo). 

Figura 9 Resultados en la etapa de pronóstico de los modelos ANFIS autónomo y MODWT-ANFIS (422). 

En relación con el segundo enfoque, el cuadro 3 muestra los resultados del cálculo de los coeficientes de correlación. En cuanto a las aproximaciones, se observa que tienen una correlación más elevada. La mayor se encuentra en el nivel uno (0.826) y van decreciendo conforme se incrementa el índice del nivel; esto indica que si el algoritmo de pronóstico puede modelar en los primeros niveles perderá menos información en el suavizado, pero la componente eliminada en altas frecuencias será mínima. Lo ideal sería que la correlación mayor existiera en un nivel más bajo (índice de nivel mayor), para que la nueva señal pudiera contener menores componentes de altas frecuencias. Por ejemplo, en el caso de niveles mensuales de lagos en Turquía, ^{Kişi (2009)} calculó el coeficiente de correlación de la componente de aproximación del nivel cinco igual a (0.890); la correlación de los detalles desciende del nivel cinco al nivel uno de 0.190 a-0.006. Por su parte, ^{Partal (2009)} trabajó con series de tiempo diarias de temperatura, radiación solar, velocidad del viento, humedad relativa y evapotranspiración en San Diego, California, EUA, y encontró un coeficiente de correlación significativo en los detalles del nivel ocho para todas esas variables. Para estos dos casos, la señal se puede reconstruir, eliminando casi por completo la componente de altas frecuencias. En la serie de la estación Ensenada será lo contrario debido a que es una serie caótica.

Cuadro 3 Coeficiente de correlación de la señal con las descomposiciones del nivel 1 a 9. 

En el cuadro 4 se presentan los resultados del modelo WN del segundo enfoque. Sólo se utiliza el algoritmo de Mallat debido a las decimaciones. El mejor resultado se logró empleando en la reconstrucción los coeficientes de aproximación y detalle del nivel seis con coeficientes de correlación de 0.305 y 0.023, respectivamente, en relación con la señal original. En otras palabras, la señal es constituida en general de altas frecuencias; por esa razón se requiere descomponerla a un nivel que corresponde a un periodo de 5.3 años y no en el nivel uno, como lo indica el coeficiente de correlación (0.826), donde la señal original y la nueva señal son más parecidas. En la figura 10 se muestran los resultados del segundo enfoque con MallatANFIS. En esa figura, los valores calculados son la línea más oscura que se observa alrededor de la media, pero en realidad es una serie de puntos que la hacen ver oscura, es decir, es la nueva señal. En esa línea, en el intervalo de pronóstico, se observa la distorsión producto del fenómeno del efecto de borde.

Cuadro 4 Resultados del modelo WN del segundo enfoque. 

Figura 10 Resultados de entrenamiento, validación y pronóstico utilizando el segundo enfoque con el algoritmo Mallat-FFBP. Se muestran los tres grupos de valores: entrenamiento, validación y pronóstico. 

Existen dos aspectos que se deben resaltar. Primero, la influencia que tiene la TW al utilizar Db4, la cual se refiere a que la nueva señal se debe a que adopta la forma de la función madre de Daubechies (comparar figuras 1a con 10), por la propiedad fractal del wavelet; es decir, es como un transgénico digital de la señal. Segundo, se genera un nuevo efecto de borde debido a que los coeficientes fueron mutilados y en el proceso de la reconstrucción esa fracción de la señal se pierde. En la figura 10 se muestran los resultados de entrenamiento, validación y pronóstico utilizando el segundo enfoque con Mallat-ANFIS. En el intervalo de pronóstico se observa la distorsión producto del fenómeno del efecto de borde.

Por otra parte, se observa en el resultado de pronóstico del cuadro 4 que utilizar FFBP para predecir la serie de tiempo en esta región es prácticamente similar a emplear la media aritmética de la señal como se esquematiza en la figura 10; además, ANFIS autónomo supera el desempeño del algoritmo Mallat-ANFIS, en términos de RMSE.

En la figura 11 se muestran los valores calculados con ANFIS autónomo y serie medida. Se observa que los valores extremos en el periodo de lluvia máxima son subestimados y en sequía se sobreestiman, aunque en ocasiones se obtienen valores negativos.

Figura 11 Resultados con el método ANFIS autónomo. Se muestran los tres grupos de valores: entrenamiento, validación y pronóstico. 

Conclusiones

En relación con los métodos autónomos estudiados, ANFIS reproduce el comportamiento de la serie de precipitación mejor que el de FFBP; esta última no logra identificar el patrón de variación de la series de precipitación; además, ANFIS es un método más sencillo que los híbridos, sin el problema del efecto de borde. En el primer enfoque, el método de pronóstico con preprocesamiento mostró mejores indicadores de desempeño en el entrenamiento, pues en el proceso de convolución, los coeficientes guardan la información futura en los instantes del periodo de entrenamiento y validación. Sin embargo, en la etapa de pronóstico influye de manera significativa el efecto de borde debido al doble pronóstico que incluye el proceso de convolución en la etapa de rellenado de la serie; esto, porque de manera involuntaria se introduce información futura, lo que es un aspecto fundamental para el pronóstico. Con las pruebas realizadas se demuestra que el coeficiente de correlación entre los datos de rellenado y los que se requieren pronosticar está directamente relacionado, y mientras más cercano a uno sea este coeficiente, mejor será el desempeño de dicho proceso; si es cercana a cero, el efecto de borde distorsiona la información y no se logra un buen pronóstico.

En el segundo enfoque se demuestra que la señal está conformada por una componente significativa de alta frecuencia; por lo tanto, los algoritmos de predicción no logran modelar hasta que la señal se aproxima a la media en los niveles más bajos y hace necesario eliminar la mayoría de las altas frecuencias en los niveles más altos. La influencia que tiene la TW al utilizar Db4, en el segundo enfoque, al encontrar la nueva señal se debe a que adopta la forma de la función madre de Daubechies. Por otra parte, se genera un nuevo efecto de borde, pues los coeficientes fueron mutilados y al momento de la reconstrucción esa fracción de la señal se pierde.

El indicador caótico de Hurst sugiere que la serie es Browniana o al azar, y no se puede predecir, lo cual es coherente con los resultados de los indicadores de desempeño en la etapa de pronóstico. Los resultados indican que la FFBP y la media de la señal muestran resultados similares (39.32 y 38.91 mm) en términos de RMSE; por otra parte, se considera que ANFIS autónomo es el más sencillo y mejor método para este tipo de series de tiempo.