Predicción de radiación solar a corto plazo mediante imágenes de cielo y visión por computador

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo de Fin de Carrera Ingeniería Electrónica, Robótica y Mecatrónica Predicción de radiación solar a corto plazo mediante imágenes de cielo y visión por computador Autor: Julián López Gordillo Tutor: Carlos Vivas Venegas Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo de Fin de Carrera Ingeniería Electrónica, Robótica y Mecatrónica Predicción de radiación solar a corto plazo mediante imágenes de cielo y visión por computador Autor: Julián López Gordillo Tutor: Carlos Vivas Venegas Profesor Contratado Doctor Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo de Fin de Carrera: Predicción de radiación solar a corto plazo mediante imágenes de cielo y visión por computador Autor: Julián López Gordillo Tutor: Carlos Vivas Venegas El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Resumen La energía solar fotovoltaica está llamada a ser una de las energías del futuro. No obstante, aun existen obstáculos por superar, como su variabilidad, especialmente aquella causada por la oclusión del sol debido a la nubosidad local. Para que pueda integrarse de forma masiva en la red eléctrica, es necesario poder predecir dichas fluctuaciones a corto plazo. Dichas predicciones se obtienen a partir de imágenes del cielo tomadas mediante cámaras a nivel terrestre, pues solo así se puede aspirar a un horizonte temporal de minutos (en comparación al horizonte temporal de horas de las imágenes por satélite). En este trabajo se realiza una revisión del estado del arte de las técnicas utilizadas en este campo de estudio, mostrando los distintos enfoques utilizados, y explicando las limitaciones de cada técnica. Para ello, se divide el problema según tres fases: segmentación semántica, estimación del movimiento, y predicción. Además, se ponen en práctica algunas de las técnicas estudiadas, mostrando cómo se integran para dar lugar a un solo algoritmo. En ese ámbito, se detalla el preprocesado realizado a las imágenes obtenidas por la cámara, además de la implementación de las demás fases ya mencionadas. Aunque la calidad de los resultados dista mucho de aquella necesaria para ser implementada en sistemas reales, los experimentos realizados constituyen una demostración didáctica de las distintas fases que forman el proceso de predicción, y pueden servir de punto de partida para aquel que busque comprender mejor dicho proceso. I Abstract So lar power is targeted towards becoming one of the main sources of energy in the future. Yet, there are still some limitations to overcome, such as its variability, especially that caused by the occlusion of the sun by the local cloudiness. For solar power to be massively integrated into the electrical grid, short-term forecasts about such fluctuations are necessary. Such forecasts are obtained from whole-sky images taken by ground-based cameras, since that is the only method able to aim at a forecast horizon of minutes (in contrast to the forecast horizon of hours obtained with satellite images). In this project, a revision of the state of the art techniques used in this field of study is conducted, showing the approaches taken, and explaining the limitations of each technique. For that matter, the process is divided in three stages: semantic segmentation, motion estimation, and forecasting. Furthermore, some of the revised techniques are put into practice, showing their integration into a single algorithm. In this case, the preprocessing applied to the images obtained by the camera is described in detail, as well as the implementation of the aforementioned stages. Although the quality of the results is far from the required to be implemented in real systems, the experiments conducted constitute a didactic demonstration of the different stages which form the forecasting process, and can serve as a starting point for those looking for understanding the process as a whole. III 2 Estado del arte En este capítulo se presentan los métodos usados en la actualidad para resolver el problema expuesto. Se supondrá cierta familiaridad con técnicas habituales en el campo de la visión por computador, por lo que no se explicará la teoría tras ellas. Se mencionarán investigaciones sobre técnicas variadas, pero se profundizará más en aquellas que han inspirado en mayor medida el presente proyecto, las cuales se pondrán en práctica en siguientes apartados. Las distintas técnicas se presentarán en un orden lógico, según la profundidad del tratamiento realizado a las imágenes. Aunque hay multitud de procedimientos y en muchos casos se prescinde de ciertas técnicas, la gran mayoría de los procedimientos siguen las tres siguientes fases: •Segmentación semántica •Estimación del movimiento •Predicción Por último, se discutirán brevemente las conclusiones que se han obtenido tras revisar y estudiar este campo. 2.1 Segmentación semántica La segmentación semántica consiste en clasificar cada píxel de la imagen según a qué objeto pertenezca. En el caso que concierne a este trabajo, las posibilidades se reducen básicamente a tres: nubes, cielo, y disco solar. La segmentación es el punto con mayor repercusión en el resultado final, y cuanto más fina sea la segmentación, menor tratamiento habrá que realizar posteriormente, y más fidedigno será el modelo del cielo con la realidad. Es por ello que también es la fase más ampliamente estudiada. Tras el proceso de segmentación, se obtiene una máscara (normalmente binaria) con la silueta de las nubes sobre el resto de objetos. De esa forma se puede extraer la información sobre qué regiones del cielo están cubiertas, y cuáles no. Es decir, esa máscara pasa a ser el modelo del cielo. 2.1.1 Selección del umbral En su forma más básica, para la obtención de dicha máscara es necesario fijar un umbral de intensidad en los píxeles de la imagen, para así separarlos de la imagen en dos conjuntos: aquellos pertenecientes a las nubes, y aquellos pertenecientes al cielo. Existen tanto algoritmos fijos (definen un único umbral consante) como adaptativos (el umbral varía según ciertas condiciones), cada uno con sus ventajas y sus desventajas. En [ 10 ], se propone un método híbrido, el HYTA (Hybrid Thresholding Algorithm, por sus siglas en inglés), que busca aunar lo mejor de ambas opciones, siendo una de las investigaciones con mayor repercusión en posteriores trabajos en el campo. Operando con los distintos canales de color, transforman la imagen de entrada en otra con un contraste más marcado 1 . Luego, categorizan dicha imagen como unimodal o bimodal, y utilizan un algoritmo 1La operación aplicada es NRBR, Normalized Red to Blue Ratio. Se explicará más sobre ella en la parte práctica. 5 6Capítulo 2. Estado del arte de segmentación fijo o de entropía cruzada mínima (MCE, por sus siglas en inglés), según interese. Los resultados experimentales validan el trabajo con una precisión en la segmentación del 88.53%, superando a cualquier otro algoritmo ideado hasta la fecha de la publicación de dicha investigación, al menos. Los mismos autores idearon un método para detectar nubes delgadas en [ 11 ]. Las nubes delgadas suponen una de las mayores dificultades en la segmentación, debido a sus bordes poco definidos y al bajo contraste con el fondo del cielo. En este caso, tratan la detección de nubes como un problema de clasificación según una red de Markov, consiguiendo unos resultados superiores a los demás algoritmos en cuanto a la detección de nubes delgadas. La idea anterior ha sido dominante en la mayoría de enfoques que la han sucedido. Sin embargo, en [ 7 ] se propone un método sin duda innovador. Las novedades introducidas en este estudio podrían resumirse como sigue: • Se comparan exhaustivamente distintos canales de color y operaciones sobre ellos, y se identifican los más adecuados para cada caso concreto. No se realiza ninguna suposición sobre los espacios de color, ni se impone ninguna restricción al tipo de imágenes de entrada. • El enfoque tomado se basa en el aprendizaje sobre un dataset, y por tanto no depende del ajuste manual de ningún parámetro o condición. Proporciona una solución sistemática. • En lugar de una clasificación binaria de las nubes, se propone una clasificación probabilística, más adecuada a la naturaleza difusa de las mismas. Ello se representa mediante máscaras de color virtual, mostrando el grado de pertenencia de un pixel concreto a la clase de nubes según su color. Aunque este trabajo promete abrir nuevas vías en la fase de segmentación, aun es pronto para ver su repercusión real. El tiempo dirá si sustituye al HYTA como enfoque dominante. 2.1.2 Redes neuronales La segmentación semántica, y el campo del tratamiento de imágenes en general, está sufriendo una revolución con la llegada de las redes neuronales, que ofrecen resultados sensiblemente superiores a las técnicas hasta ahora utilizadas, incluidas muchas del aprendizaje automático clásico. En [ 20 ] se realiza una comparación entre un método clásico de aprendizaje automático, las máquinas de vectores de soporte (SVM, Support Vector Machine), y una red neuronal. Mientras que en ciertos casos ambos métodos proporcionaban resultados similares, en aquellos donde la SVM fallaba rotundamente, la red neuronal conseguía funcionar con normalidad, proporcionando unos resultados consistentes. En este otro estudio, [ 22 ], la detección de nubes se enfoca como un problema de clasificación mediante una red neuronal de tres capas. Como entrada, a la red le fueron proporcionadas los canales RGB de las imágenes del cielo por separado, pero también canales adicionales como el resultante de la ya mencionada imagen NRBR. Se etiquetaron manualmente 600 imágenes como parte del dataset de aprendizaje de la red, la cual produce como salida una clasificación de cada píxel como parte del cielo, del disco solar, o de las nubes. Los resultados muestran que el uso de redes neuronales permite obtener unos algoritmos de segmentación más inteligentes que los tradicionales basados en umbrales, pues los primeros son capaces de evitar errores de mayor sutileza como los reflejos de lente, entre otros. 2.2 Estimación del movimiento Una vez las imágenes son obtenidas y procesadas para obtener la máscara de nubes, dichas imágenes se consideran en una secuencia para poder inferir el movimiento dominante de las masas nubosas, y así disponer de información sobre el futuro con la que realizar una predicción. El método de estimación del movimiento tiene que ir en consonancia al modelo del cielo empleado, es decir, a cómo se presenta la información en las imágenes tras la segmentación. Por ejemplo, en ciertos casos se realiza una transformación geométrica para corregir la distorsión introducida por las lentes de gran angular, mientras que en otros se trabaja directamente con las imágenes obtenidas. En la estimación del movimiento, dichas distorsiones deben ser consideradas, y ciertos métodos se comportan mejor que otros dependiendo de ese y otros muchos factores, como se verá a continuación. 2.2 Estimación del movimiento 7 2.2.1 Block matching Las técnicas de block matching se basan en dividir la imagen en un conjunto de bloques, y encontrar la mejor correspondencia entre los bloques de una imagen y los de la siguiente imagen de la secuencia según un criterio de semejanza. Al utilizar secciones locales de cada imagen, se trata de un método robusto frente al ruido y a variaciones de brillo dentro de la propia imagen. Puesto que el movimiento de las nubes es mayoritariamente de traslación (sin cambios de escala, aunque sí pueden existir deformaciones), este método puede representar verídicamente el movimiento de las nubes, como se demostró en [9]. Sin embargo, la mayoría de los algoritmos de block matching emplean bloques de tamaño y posición predefinidos, en forma de una malla en la imagen, [ 5 ]. Dado que las nubes pueden tomar muy variadas formas, tamaños, y posiciones en la imagen, el rendimiento de este tipo de métodos no es consistente y varía dependiendo de lo que se ajusten las características anteriores a los bloques escogidos. Para solventar dichos problemas se han propuesto algunas mejoras, como la asignación dinámica de bloques según la máscara de las nubes en la imagen del cielo, [ 17 ], o una estructura jerárquica de bloques para estimar el movimiento con mayor precisión. Aunque dichas modificaciones consiguen mejorar los resultados, la raíz del problema permanece, pues la base del método es una división poco flexible de la imagen por bloques cuya forma es invariable. 2.2.2 Flujo óptico En el campo de la visión por computador, el problema de la estimación de movimiento suele resolverse mediante flujo óptico. Los algoritmos de flujo óptico atacan el problema desde un nivel más bajo que los de block matching. Estos algoritmos consiguen estimar el movimiento a nivel de píxel, y para ello buscan minimizar una función objetivo (usualmente una función en términos de energía), asumiendo continuidad en el brillo entre imágenes. En comparación con los métodos de block matching, estos resultan significativamente más flexibles. A cambio, son más sensibles al ruido y a variaciones en el brillo, ambas características comunes en las imágenes del cielo. De entre las diferentes implementaciones de flujo óptico, la más utilizada en este campo es la de Lucas- Kanade, [13]. Para atenuar la influencia del ruido, en [22] se propone un filtro de similitud entre vecindades del campo de velocidades. Se asume un mínimo de similitud entre vectores de una misma vecindad, y luego se insiste en dicha similitud entre vecindades colindantes. Con ello se consiguen disminuir las direcciones erróneas en el campo vectorial, así como unos resultados más consistentes. La deformación de las nubes y los cambios en el brillo suponen un reto para este algoritmo. En [ 23 ] se toma la ya mencionada imagen NRBR como entrada para el algoritmo de Lucas-Kanade, consiguiendo así unas características más estables para ser rastreadas por el algoritmo en sucesivas iteraciones, pues las variaciones en el brillo quedan suavizadas. Ambas ideas serán útiles en la parte práctica del presente trabajo. 2.2.3 Otros métodos Aunque las técnicas de block matching y flujo óptico cubren la mayoría de los métodos utilizados para estimar el movimiento, existen otros que resuelven el problema de una forma distinta. Por ejemplo, en [ 1 , 2 ], se propone un registro de imágenes no rígido para poder abordar la dinámica del movimiento de las nubes y su deformación en el tiempo. De esta forma, consiguen una gran efectividad en aquellos aspectos en los que las técnicas de block matching son más débiles. Para ello se basan en el algoritmo "Demons" ideado por Thyrion, [21], el cual es ampliamente utilizado en imágenes médicas. 8Capítulo 2. Estado del arte 2.3 Predicción Por último, tras la estimación del movimiento viene la fase de predicción. Existen multitud de técnicas que, en base a los datos obtenidos en fases anteriores (y posiblemente fusionando otros datos nuevos), realizan una predicción sobre la posición de las nubes, el estado del disco solar, o la irradiancia recibida. Sin duda, esta fase es donde más diversidad de técnicas se han ideado. Desde técnicas tradicionales de aprendizaje automático como las ya mencionadas SVM, a métodos diseñados específicamente para este tipo de casos. Aunque en los últimos años las redes neuronales están dominando en el ámbito de la predicción, cabe destacar el siguiente [23] por el ingenioso uso de otras técnicas ya mencionadas: En él, usan redes neuronales para la detección de nubes, y Lucas-Kanade para la estimación del movimiento. Utilizando el campo de velocidades obtenido, junto con las regiones clasificadas como nubes, se genera una predicción sobre qué nubes cruzarán su trayectoria con el sol. Dicha predicción se presenta en la forma de un tren de ondas, para poder ser visualizada en el tiempo de una forma clara. Además, en la fase de predicción es habitual incorporar datos adicionales que puedan complementar la información obtenida por las cámaras. Típicamente, se suele utilizar la medida de Irradiancia Global Horizontal (GHI, por sus siglas en inglés) o Irradiancia Normal Directa (DNI). En [ 19 ], por ejemplo, se combina una red neuronal de tres capas y una fase de preprocesado basada en lógica borrosa para realizar predicciones en base a los datos de irradiación obtenidos, combinando así distintas fuentes de información. 2.4 Conclusiones Tras estudiar una gran variedad de técnicas en las distintas fases del proceso, se pueden obtener ciertas conclusiones sobre los retos que aun quedan por superar, así sobre los enfoques más acertados para enfrentarse a ellos. Una de las mayores dificultades en este campo, sin duda alguna, reside en lo cambiante de las condiciones a las que un sistema debe trabajar. El sol se desplaza a lo largo del cielo, pudiendo pasar de iluminar directamente a la cámara a no hacerlo, con los cambios de intensidad luminosa que ello conlleva. Además, la trayectoria aparente del sol varía a lo largo del año. El clima local también afecta, por supuesto, y ciertos días estará más nublado que otros. Y por supuesto, el momento en el que el sol es obstruido por las nubes, las condiciones lumínicas cambian drásticamente, lo cual ha de ser tenido en cuenta en la calibración de parámetros. Además, los tamaños, formas, velocidades y direcciones de las nubes también son altamente variables, llegando incluso a coexistir más de un tipo de nube en el cielo al mismo tiempo, en distintas capas. Cada capa posee una velocidad independiente, y dado que la mayoría de algoritmos no están preparados para considerar más de una capa, este hecho puede falsear por completo cualquier medida obtenida. Todo ello dificulta el correcto funcionamiento de ciertos algoritmos, que a menudo parten de unas hipótesis de constancia en ciertas características (por ejemplo, constancia en el brillo en el algoritmo de flujo óptico de Lucas-Kanade). Cada técnica tiene un escenario ideal en el que su rendimiento es máximo, pero también otras condiciones en las que sus resultados se deterioran sensiblemente. Es por ello que se tiende a la hibridación de algoritmos, con el objetivo de obtener las ventajas de unos y otros, y poder así lidiar con situaciones cambiantes. Muchas de las técnicas mencionadas hacen uso de métodos híbridos y adaptativos para poder funcionar con precisión en un amplio abanico de situaciones. No obstante, las redes neuronales han demostrado ser notablemente más flexibles que los algoritmos tradicionales, y precisamente su campo de excelencia es el del reconocimiento de patrones, particularmente en imágenes. Todo apunta a que dichas técnicas serán las que dominarán este campo en el futuro, aunque a corto plazo lo más probable es que sean emparejadas con otras técnicas de forma complementaria, como hemos podido observar en los estudios considerados en este capítulo. 3 Trabajo práctico En este capítulo pondremos en práctica ciertas técnicas de las presentadas en el capítulo anterior. En primer lugar, se expondrá un ejemplo de las imágenes con las que se ha trabajado, y a continuación, paso por paso, se detallará el tratamiento realizado a dichas imágenes. Luego, se explicará cómo se ha obtenido información a partir de ellas, y por último, se someterá al algoritmo a una batería de pruebas para evaluar su comportamiento. Las distintas fases del proceso se estructuran de la siguiente forma: •Obtención de imágenes •Segmentación •Operaciones morfológicas •Estimación del movimiento •Predicción 3.1 Obtención de imágenes La obtención de las imágenes se realiza mediante una cámara situada en lo alto de uno de los edificios de los laboratorios de la Escuela Técnica Superior de Ingeniería (ETSI), enfocando al cielo. Un ejemplo de imagen obtenida con dicha cámara es el de la Figura 3.1. La resolución de las imágenes es de 2048x1536 píxeles. Como se explicó en el capítulo anterior, en ocasiones se opta por corregir el efecto producido por la lente de gran angular mediante una transformación geométrica. En este caso, no se ha realizado tal procedimiento, y se ha trabajado directamente sobre las imágenes con dicho efecto. 9 10 Capítulo 3. Trabajo práctico Figura 3.1 Ejemplo de imagen obtenida por la cámara. 3.2 Segmentación semántica 11 3.2 Segmentación semántica La primera fase del proceso es la de segmentación semántica. El modelo del cielo consistirá en una máscara binaria, cuyos píxeles de baja intensidad (negro) serán los correspondientes a las nubes, y los de alta (blanco), al cielo. La segmentación podría intentar realizarse directamente sobre la imagen obtenida por la cámara, pero aunque a la vista humana le resulte sencillo, al descomponer la imagen en los tres canales (rojo, verde y azul) la diferencia no resulta tan evidente. En numerosas investigaciones, como la ya mencionada [ 10 ], se muestra que la segmentación resulta más evidente al usar la diferencia normalizada de los niveles de rojo y azul 1 , NRBR (Normalized Red-Blue Ratio), la cual se define como: NRBR = B−R B+R Al descomponer la imagen original en sus componentes verde, azul, y roja, y realizar la operación indicada (donde NRBR es la imagen resultante, Bhace referencia al canal azul y Ral rojo), se obtiene la imagen de la Figura 3.2. Figura 3.2 Imagen obtenida tras aplicar la operación NRBR. En ella se puede observar un mayor contraste entre las masas nubosas y el cielo. 3.2.1 Método de Otsu Una vez obtenida la imagen NRBR, para la obtención de un umbral con el que realizar la segmentación, se aplicó el criterio de Otsu como primera aproximación. Dicho método dio como resultado un umbral θ=0.0667 (en la escala de intensidad de 0-255, se correspondería con un valor de 17, aproximadamente), con el que se obtuvo la máscara binaria que se muestra en la Figura 3.3. Debido a que los valores de intensidad correspondientes a las nubes no son uniformes a lo largo de la imagen, es imposible segmentar a la perfección mediante un sólo umbral. Por tanto, es necesario llegar a una solución de compromiso entre los píxeles clasificados erróneamente como parte del cielo y como parte 1 Esto se debe a la llamada Dispersión de Rayleigh, que describe la dispersión de la radiación electromagnética por partículas cuyo tamaño es mucho menor que la longitud de onda de los fotones dispersados. Dicho efecto es el principal causante de que el cielo se vea azul. 12 Capítulo 3. Trabajo práctico Figura 3.3 Máscara binaria según el criterio de Otsu. de las nubes, y corregir en una fase posterior dichas clasificaciones erróneas. Al experimentar con distintos umbrales, se consiguieron resultados más precisos al proporcionado por el criterio de Otsu. La solución más favorable se encontró para el valor θ=0.085 ( 21 en la escala 0-255), en la Figura 3.4 se muestra la máscara resultante. Figura 3.4 Máscara binaria obtenida empíricamente. 3.2.2 Exclusión del disco solar Nótese que el disco solar consta de unos niveles de intensidad muy elevados, al igual que en las regiones pertenecientes a las nubes. Dado que no es posible diferenciar ambos elementos según su intensidad, ha de hacerse por algún otro método. En este caso, se ha optado por excluir el disco solar de la máscara de nubes, 3.2 Segmentación semántica 13 para obtener así una máscara más precisa. La exclusión se ha realizado manualmente, seleccionando un área lo suficientemente grande como para que contenga al disco solar, y también a gran parte de los destellos causados en la lente, como puede observarse en la Figura 3.5. El resto de los destellos serán corregidos posteriormente. Figura 3.5 Máscara binaria con el disco solar excluido. 20 Capítulo 3. Trabajo práctico 3.4.2 Block-matching En el caso del block matching, en primer lugar hay que seleccionar un tamaño apropiado para los bloques, y para el máximo desplazamiento aplicado a la hora de buscar en la siguiente imagen de la secuencia. Dichos parámetros dependerán, entre otras cosas, de la resolución de la imagen, del tamaño de los objetos a rastrear, de las restricciones temporales sobre el tiempo de computación, etc. Dada la resolución de la imagen, y el tamaño medio de las nubes, se ha considerado que 50 píxeles es un valor razonable para ambos parámetros. Los resultados obtenidos se muestran en la Figura 3.14. Figura 3.14 Campo de vectores producido por el algoritmo de Block-matching. Se puede observar que aunque la presentación de los vectores es distinta con respecto al método de flujo óptico, hay coherencia entre ambos métodos en cuanto al vector global de movimiento. Con ello en mente, se ha optado por utilizar el método de flujo óptico, al ser la información sobre el movimiento más intuitiva en ese caso. 3.5 Predicción 21 3.5 Predicción El objetivo fundamental del proceso es obtener una predicción sobre el estado del cielo, y particularmente del disco solar. Para ello, una vez se ha obtenido el vector global de movimiento, se pueden trasladar toda la masa nubosa según la dirección del vector, y en una magnitud acorde al horizonte temporal en el que se quiera realizar la predicción. 3.5.1 Elección del horizonte temporal Como se comentó en la introducción, una de las ventajas de realizar este proceso mediante cámaras situadas a nivel terrestre y no mediante satélites es la posibilidad de obtener predicciones de corto plazo, del orden de minutos. Con eso en mente, el horizonte temporal máximo es aquel que transcurre entre que las nubes objetivo acaban de entrar en la imagen y hasta que llegan al disco solar. En las condiciones de este proyecto, se ha comprobado que dicho horizonte máximo oscila entre los 4 y 8 minutos, dependiendo de la posición del sol (hora del día), y la dirección y velocidad de las nubes respecto del mismo. Puesto que el objetivo de esta sección es el de ilustrar el proceso en su totalidad, y tras experimentar con distintos intervalos de tiempo, el horizonte temporal escogido para esta sección es de 3 minutos. Dicho horizonte suele proporcionar unos resultados equilibrados en cuanto a la calidad de la predicción y su lejanía en el horizonte temporal (Tabla 3.1). No obstante, en la siguiente sección se mostrará un conjunto más amplio de resultados donde el horizonte temporal se supone variable, para así poder observar cómo se modifica la precisión de las predicciones realizadas conforme el horizonte temporal aumenta. 3.5.2 Obtención de resultados Una vez realizada la traslación de las nubes, puede calcularse el porcentaje de cielo cubierto por nubes como la proporción de píxeles clasificados como nubes sobre el total en la zona central de la imagen. Del mismo modo, se puede comprobar la presencia de nubes en la región del disco solar, para poder determinar su estado. Se han considerado tres estados posibles: despejado, cubierto, y parcialmente cubierto. En la Figura 3.15 se puede observar la predicción realizada para la imagen, con un horizonte temporal de 3 minutos. En la parte superior de la imagen se encuentra la predicción sobre el estado y el porcentaje del cielo cubierto. Figura 3.15 Predicción realizada mostrando énfasis en la región del disco solar. 22 Capítulo 3. Trabajo práctico Pueden distinguirse dos círculos en la imagen: uno abarca ampliamente la región central de la imagen, y el otro, de menor tamaño, está situado aproximadamente en el centro de la imagen en este caso. El primero constituye la mayor parte del cielo observable (las regiones periféricas contienen a los edificios y están demasiado distorsionadas, por lo que no se consideran), y es la región sobre la que se calcula el porcentaje de cielo cubierto. El otro círculo, de menor tamaño y situado cerca del centro de la imagen, representa al disco solar. Sirve de referencia para evaluar si dicha región está ocluida por nubes o no. En este caso, puede observarse que el disco solar está mayoritariamente despejado, como así lo indica la predicción, y puede verificarse cualitativamente la predicción realizada sobre el porcentaje de cielo cubierto. Nótese que al realizar la predicción sólo se traslada la región central de la imagen, y que al hacerlo se originan incoherencias en la región periférica. No obstante, las regiones periféricas no son de nuestro interés, y las imperfecciones en dichas regiones no afectan sensiblemente a nuestro objetivo. Lo que sí ha de tenerse en cuenta es que, al trasladar dicha región central de la imagen, queda una zona "vacía" en ciertas partes que antes ocupaba la imagen. Dichas zonas se rellenan con nubes ante la imposibilidad de saber si en dichas zonas entrarán o no nubes en el futuro (esa información queda fuera de nuestro horizonte temporal). Este efecto supone una fuente de error, y es mayor cuando mayor sea el horizonte temporal escogido. En este punto, ya se ha obtenido información sobre el estado futuro del disco solar, el dato más relevante en lo que concierne a los paneles fotovoltaicos, por ejemplo. Además, también se ha obtenido una medida del porcentaje de cielo cubierto, información complementaria que también podría ser tenida en cuenta a la hora de tomar decisiones. Por tanto, la predicción queda realizada. 3.6 Resultados 23 3.6 Resultados A continuación se mostrarán varios experimentos llevados a cabo adicionalmente, para la comprobación del funcionamiento del algoritmo propuesto. Para ello, se dispone de una batería de imágenes del cielo tomadas durante varios días consecutivos. Esto permite comparar las predicciones realizadas sobre el estado del cielo en un momento dado con el estado del cielo real que se registró en ese momento, pudiendo así obtener una métrica de la precisión de nuestro algoritmo. En primer lugar se mostrarán una serie de experimentos ilustrativos. Dichos experimentos han sido seleccionados para mostrar tanto aquellas situaciones en las que el algoritmo se comporta adecuadamente como aquellas en las que falla, y explicar las razones tras dicho funcionamiento. Luego, se mostrarán y comentarán los resultados obtenidos tras someter al algoritmo a una batería intensiva de pruebas, en las que los tipos de nubes, la posición del sol, y el horizonte temporal son variables. Los resultados obtenidos serán representados en forma de gráficas para una interpretación más intuitiva. 3.6.1 Experimentos ilustrativos Se han seleccionado varias secuencias de imágenes con distintos tipos de situaciones y nubes. En todos ellos, las nubes están en camino de cubrir al disco solar, pero el tamaño y forma de las nubes varía, y con ellos, la calidad de los resultados. Se mostrarán tres casos típicos de nubes en orden creciente de complejidad, y cómo el algoritmo se desenvuelve en dichos escenarios. Caso favorable En primer lugar, se ha elegido una imagen en la que se aproxima una masa nubosa de gran tamaño y bien definida (Figura 3.16). Figura 3.16 Imagen del primer caso seleccionado. Nubes de este tipo suponen los casos más sencillos ya que ocluyen completamente al disco solar, por lo que se puede esperar una predicción certera sobre su estado. Además, las predicciones sobre el porcentaje 24 Capítulo 3. Trabajo práctico cubierto del cielo resultan razonables en general, como puede observarse en las siguientes páginas, en la Figura 3.17. Como medida de la precisión de la predicción realizada, se presenta la imagen real tomada 3 minutos después (Figura 3.18). Caso problemático El siguiente caso es un ejemplo claro de las limitaciones del método utilizado. Se trata de un conjunto de numerosas nubes de pequeño tamaño y con poca separación entre ellas (Figura 3.21). Aunque se encuentran bien definidas, el algoritmo no realiza un buen trabajo debido a las operaciones morfológicas realizadas, que cubren muchos de los huecos entre las nubes (Figura 3.20). Se genera una predicción imprecisa sobre el estado del disco solar, aunque al menos estará por el lado de la seguridad. Caso inalcanzable El último caso es en general el más desafiante para la mayoría de algoritmos revisados en la literatura. Se trata de nubes finas y translúcidas, vagamente definidas. En la Figura 3.21 resulta evidente que dichas nubes ni siquiera consiguen ocluir por completo los rayos solares al cubrir el sol, y sin embargo engañan a la mayoría de algoritmos en cuanto a su detección. En este caso se espera un fracaso total del algoritmo, al ser incapaz de clasificar con precisión las delgadas masas nubosas, y mucho menos de saber que estas no impedirán la penetración de los rayos solares al cubrir al sol. El resultado se muestra en la Figura 3.22. 3.6 Resultados 25 Figura 3.17 Predicción realizada sobre el primer caso seleccionado, con un horizonte temporal de 3 minutos. Figura 3.18 Imagen real (tras el procesado) a los 3 minutos de la anterior. 26 Capítulo 3. Trabajo práctico Figura 3.19 Imagen del segundo caso seleccionado. Figura 3.20 Predicción realizada sobre el segundo caso seleccionado, con un horizonte temporal de 3 minutos. 3.6 Resultados 27 Figura 3.21 Imagen del tercer caso seleccionado. Figura 3.22 Predicción realizada sobre el tercer caso seleccionado, con un horizonte temporal de 3 minutos. 28 Capítulo 3. Trabajo práctico 3.6.2 Baterías de pruebas En los casos anteriores, se seleccionaron casos específicos, con intervalos cortos de tiempo, suficientes para mostrar los resultados concretos que se buscaban explicar. A continuación, sin embargo, se tomarán varias baterías de pruebas de larga duración, y se ejecutará el algoritmo sobre ellas. Se utilizarán un total de 4 baterías de pruebas, las cuales consisten en un conjunto de imágenes del cielo tomadas cada 5 segundos, durante 1 hora, en distintos momentos del día (unas 720 imágenes por batería). Tal lapso de tiempo es suficiente para que la nubosidad percibida de forma local varíe sensiblemente, pudiendo tener en cada batería numerosos conjuntos de nubes que entran y salen de la imagen, incluso transiciones en el estado del cielo (de totalmente despejado a cubierto, o viceversa). Cada batería se ejecuta un total de 6 veces, variando el horizonte temporal de 1 a 6 minutos. Ello supone un total de 24 ejecuciones sobre las baterías de pruebas, en las cuales, como se ha argumentado, tanto las condiciones de nubosidad como los horizontes temporales son variables. Por tanto, los resultados obtenidos deberían de constituir una muestra aceptablemente fiable del comportamiento general del algoritmo, y deberían de arrojar luz sobre la influencia del horizonte temporal en el proceso de predicción. Gráficas Dicho esto, en las siguientes páginas se muestran los resultados obtenidos. Para cada ejecución, se ha calculado el error cuadrático medio en las predicciones realizadas durante toda la ejecución. Luego, para cada batería, se han agrupado dichos errores en gráficas, mostrando la evolución conforme se aumenta el horizonte temporal. En líneas generales, los valores y tendencias del error cometido en las distintas baterías de pruebas se mantiene. No obstante, como se puede observar en la Figura 3.24, en la segunda batería de pruebas el error cometido es considerablemente más alto. Ello se debe a la presencia de nubes delgadas, que, como se explicó anteriormente, suponen un obstáculo insalvable para el algoritmo desarrollado. En cuanto a la evolución de los resultados con el horizonte temporal, el algoritmo pierde precisión conforme aumenta el mismo, como es de esperar. Hasta los 3 minutos, el error es admisible, por lo general, y la predicción es lo suficientemente lejana en el tiempo como para que se produzcan cambios en el cielo. Supone el valor de compromiso entre el horizonte en el que el algoritmo se desenvuelve bien, y un intervalo de tiempo lo suficientemente extendido como para que la predicción sea relevante, y haya una ventana de tiempo útil en la que se pueda actuar. 3.6 Resultados 29 Figura 3.23 Resultados en las predicciones sobre la batería de pruebas nº1. Figura 3.24 Resultados en las predicciones sobre la batería de pruebas nº2. Apéndice A Código desarrollado En este apéndice se recoge el código desarrollado en el presente proyecto, de forma que el trabajo pueda ser continuado desde donde se concluye aquí sin demasiados obstáculos. A.1 Software utilizado El proyecto ha sido realizado en MATLAB®R2018a, y se ha hecho uso de los siguientes toolbox: •Computer Vision System Toolbox •Image Processing Toolbox •Curve Fitting Toolbox En las siguientes secciones se muestra el código comentado, acompañado de explicaciones sobre cómo ejecutarlo. A.2 Script principal El siguiente código recoge todas las funcionalidades realizadas por el algoritmo, a excepción del preprocesado, que se tratará en la siguiente sección. Para una correcta ejecución del código, habría que introducir la ubicación de la base de datos de imágenes, así como cargar las máscaras que representan al disco solar y al cielo relevante en el archivo ’TFG_ vars.mat’. Además, según el horizonte temporal en el que se quiera realizar la predicción, y el intervalo de toma de imágenes, han de modificarse las variables saltos yfactor_ desp, las cuales influyen en el proceso de predicción, y han de ser calibradas experimentalmente. Una vez solventados dichos puntos, el algoritmo está preparado para funcionar automáticamente. Sin más dilación, a continuación se muestra el código: Código A.1 Código ejecutable para el proceso de predicción. % Código del script principal %% Inicialización de variables y datos clearvars; clc; close all; % Carga de parámetros en el workspace. Se incluyen las máscaras binarias % seleccionando el disco solar, las cuales han sido creadas manualmente y % dependen de las imágenes usadas. Dicho archivo se usa como base de datos % para los resultados, también. load(’tfg_vars.mat’); 37 38 Capítulo A. Código desarrollado % Creación de la base de datos de las imágenes (con función de % preprocesado) skyds = imageDatastore(’D4-1430_1530’,’ReadFcn’,@DatastorePreprocessing); % Objeto de flujo óptico (método de Lucas-Kanade) opticFlow = opticalFlowLK(’NoiseThreshold’,0.09); % Lectura de imágenes y carga en el workspace numFrames = numel(skyds.Files); I = read(skyds); sequence = zeros([size(I) numFrames],class(I)); sequence(:,:,1) = I; for p = 2:numFrames I = read(skyds); sequence(:,:,p) = I; end % No ¯de saltos entre la imagen actual, y aquella que se intenta % predecir (depende del intervalo de captura de imágenes y del % horizonte temporal de predicción) saltos = 72; % Factor de desplazamiento para la traslación en el proceso de % predicción (varía según el horizonte temporal) factor_desp = 64; % Variable auxiliar, contador de predicciones t = 0; % Oclusión real del cielo occlusion_sky_r = zeros(ceil((numFrames-1)/3),1); % Oclusión predicha del cielo occlusion_sky_p = zeros(ceil((numFrames-1)/3),1); % Variable auxiliar para el cómputo del vector global de movimiento global_motion_aux = zeros(ceil((numFrames-1)/3),2); % Error cuadrático medio cometido err = zeros(ceil((numFrames-1)/3),2); %% Flujo óptico % Estimación del flujo óptico en la secuencia figure(1) for k = 1:3:numFrames if k == 1 imshow(sequence(:,:,k)); title(sprintf(’Imagen inicial’)); flow = estimateFlow(opticFlow,sequence(:,:,k)); pause(1); elseif k < numFrames-saltos t = t+1; flow = estimateFlow(opticFlow,sequence(:,:,k)); % Corrección de ruido % Vectores del campo vectorial de flujo óptico anómalos son % eliminados, pues se deben a la formación/evaporación de nubes Vxf = flow.Vx; Vyf = flow.Vy; A.2 Script principal 39 [M,N] = size(Vxf); Magf = zeros(M,N); Vxt = zeros(M,N); Vyt = zeros(M,N); Vx_sum = double(0); Vy_sum = double(0); cont = 0; for i=1:M for j=1:N Magf(i,j) = sqrt(Vxf(i,j)^2 + Vyf(i,j)^2); if Magf(i,j) > 2 Vxf(i,j) = 0; Vyf(i,j) = 0; end % Cálculo del vector global de movimiento if Magf(i,j) > 0.1 && Magf(i,j) < 2 Vx_sum = Vx_sum + double(Vxf(i,j)); Vy_sum = Vy_sum + double(Vyf(i,j)); cont = cont + 1; end end end Vxt(831,1056) = (Vx_sum/cont)*50; Vyt(831,1056) = (Vy_sum/cont)*50; % Suavizado del vector global de movimiento. Se promedian los % vectores globales de movimiento entre una vecindad de imágenes % para reducir el impacto de cambios bruscos en el brillo. Este % suavizado no desvirtúa los datos, pues la velocidad de las nubes % no cambia de forma brusca (mientras que el brillo sí, debido a la % oclusión del sol, distorsionando el flujo óptico) if t == 1 global_motion_aux(t,:) = [Vxt(831,1056) Vyt(831,1056)]; global_motion_vtr = [Vxt(831,1056) Vyt(831,1056)]; else global_motion_aux(1:t-1,:) = global_motion_aux(2:t,:); global_motion_aux(t,:) = [Vxt(831,1056) Vyt(831,1056)]; global_motion_vtr = mean(global_motion_aux(1:t,:)); end Vxt(831,1056) = global_motion_vtr(1); Vyt(831,1056) = global_motion_vtr(2); % Descomentar para visualizar el flujo óptico corregido % imshow(sequence(:,:,k)); % title(sprintf(’Imagen #%d’,k)); % flow3 = opticalFlow(Vxt,Vyt); % flow2 = opticalFlow(Vxf, Vyf); % hold on; % plot(flow2,’DecimationFactor’,[5 5],’ScaleFactor’,24); % plot(flow3,’DecimationFactor’,[5 5],’ScaleFactor’,24); % hold off; %% Predicción % Imágenes actual y futura fi = sequence(:,:,k); ff = sequence(:,:,k+saltos); % Traslación de la región central de la imagen según el horizonte 40 Capítulo A. Código desarrollado % temporal escogido. ’FillValues’ rellena las regiones de % incertidumbre con nubes (0). filt = @(binary_image) imtranslate(binary_image,factor_desp* global_motion_vtr,... ’FillValues’,0); gf = roifilt2(fi, mask_sky, filt); % Métricas % Cálculo de la oclusión del disco solar stats = regionprops(mask_sun,gf,’Area’,’MeanIntensity’); oc_correction = 255 - stats.MeanIntensity; occlusion_sun = oc_correction*100/255; if occlusion_sun < 20 estado_p = ’\color{blue}Despejado’; elseif occlusion_sun > 80 estado_p = ’\color{red}Cubierto’; else estado_p = ’\color{orange}Parcialmente cubierto’; end % Oclusión del cielo predicha stats = regionprops(mask_sky,gf,’Area’,’MeanIntensity’); oc_correction = 255 - stats.MeanIntensity; occlusion_sky_p(t) = oc_correction*100/255; % Oclusión del cielo real stats = regionprops(mask_sky,ff,’Area’,’MeanIntensity’); oc_correction = 255 - stats.MeanIntensity; occlusion_sky_r(t) = oc_correction*100/255; % Para la visualiación de las predicciones realizadas y la imagen % real, descomentar las siguientes líneas % figure(2); imshowpair(gf,mask_ss,’blend’) % title({[’Estado del disco solar: ’,estado_p, ... % ’. Porcentaje de cielo cubierto: ’,num2str(occlusion_sky_p(t)), ’%’]},’FontSize’,18); % figure(3); imshowpair(ff,mask_sky,’blend’) % title([’Oclusión solar real: ’,num2str(occlusion_sky_r), ’%. Estado: ’,estado_r]) % title({[’Porcentaje de cielo cubierto real: ’,num2str(occlusion_sky_r) , ’%’]},’FontSize’,18) % Descomentar en caso de que se realice alguna visualización de % imágenes % pause(1); % Cálculo del error cuadrático medio cometido err(t) = immse(occlusion_sky_p(t), occlusion_sky_r(t)); end end %% Presentación y conservación de los datos obtenidos % Eliminación de elementos sobrantes del vector A.3 Preprocesado de imágenes 41 err = nonzeros(err’); mean(err) % Los datos de interés pueden ser guardados para su análisis o futuras % pruebas adicionales save(’tfg_vars.mat’,’global_motion_vtr’,’err’,’-append’); A.3 Preprocesado de imágenes Previamente a su carga en la base de datos de imágenes y en el workspace, las imágenes son preprocesadas, y ello se realiza mediante una función auxiliar que se invoca en la creación del datastore. El código de dicha función se muestra a continuación: Código A.2 Función de preprocesado de las imágenes de la base de datos. function go = DatastorePreprocessing(filename) % Función de preprocesado de imágenes almacenadas en el datastore % Carga la máscara del disco solar al workspace. Dichas máscaras han sido % realizadas manualmente sobre las imágenes load(’tfg_vars’,’mask_sun’); f = imread(filename); [M, N, ~] = size(f); % Obtención de los canales de rojo y azul de la imagen fr = double(f(:,:,1)); fb = double(f(:,:,3)); % Operación NRBR (Normalised Red to Blue Ratio) f_NRBR = (fb-fr)./(fb+fr); % Adecuación del resultado al formato de imagen for i=1:M for j=1:N if f_NRBR(i,j)<0 || isnan(f_NRBR(i,j)) f_NRBR(i,j) = 0; end end end f_NRBR = 255*f_NRBR; g = uint8(f_NRBR); %% Creación de la máscara binaria thr = 0.085; gb = or(im2bw(g, thr), mask_sun); % Exclusión del disco solar gb = uint8(255*gb); go = gb; %% Operaciones morfológicas % Creación de los elementos estructuradores SEe = strel(’disk’,7); 42 Capítulo A. Código desarrollado SEc = strel(’disk’, 20); % Erosión y cierre gbe = imerode(gb,SEe); gbc = imclose(gbe,SEc); go = gbc; Índice de Figuras 3.1 Ejemplo de imagen obtenida por la cámara 10 3.2 Imagen obtenida tras aplicar la operación NRBR 11 3.3 Máscara binaria según el criterio de Otsu 12 3.4 Máscara binaria obtenida empíricamente 12 3.5 Máscara binaria con el disco solar excluido 13 3.6 Erosión con elemento estructurador en forma de diamante 14 3.7 Erosión con elemento estructurador en forma de disco circular 15 3.8 Erosión con elemento estructurador en forma de disco circular con tamaño reducido 15 3.9 Contorno añadido por la operación de erosión (verde) 16 3.10 Contorno eliminado por la operación de cierre (magenta) 16 3.11 Máscara binaria definitiva 17 3.12 Resultados del algoritmo de flujo óptico de Lucas-Kanade 18 3.13 Campo de vectores filtrado producto de Lucas-Kanade (azul), y vector global de movimiento (naranja) 19 3.14 Campo de vectores producido por el algoritmo de Block-matching 20 3.15 Predicción realizada mostrando énfasis en la región del disco solar 21 3.16 Imagen del primer caso seleccionado 23 3.17 Predicción realizada sobre el primer caso seleccionado, con un horizonte temporal de 3 minutos 25 3.18 Imagen real (tras el procesado) a los 3 minutos de la anterior 25 3.19 Imagen del segundo caso seleccionado 26 3.20 Predicción realizada sobre el segundo caso seleccionado, con un horizonte temporal de 3 minutos 26 3.21 Imagen del tercer caso seleccionado 27 3.22 Predicción realizada sobre el tercer caso seleccionado, con un horizonte temporal de 3 minutos 27 3.23 Resultados en las predicciones sobre la batería de pruebas nº 1 29 3.24 Resultados en las predicciones sobre la batería de pruebas nº 2 29 3.25 Resultados en las predicciones sobre la batería de pruebas nº 3 30 3.26 Resultados en las predicciones sobre la batería de pruebas nº 4 30 43 Bibliografía [1] David Bernecker, Christian Riess, Elli Angelopoulou, and Joachim Hornegger, Towards Improving Solar Irradiance Forecasts with Methods from Computer Vision, Computer Vision in Applications Workshop (DAGM), 2012. [2] , Continuous short-term irradiance forecasts using sky images, Solar Energy (2014). [3] ,Continuous short-term irradiance forecasts using sky images, Solar Energy 110 (2014), 303– 315. [4] Viv Bone, John Pidgeon, Michael Kearney, and Ananthanarayanan Veeraragavan, Intra-hour direct normal irradiance forecasting through adaptive clear-sky modelling and cloud tracking, Solar Energy 159 (2018), no. October 2017, 852–867. [5] Chi Wai Chow, Bryan Urquhart, Matthew Lave, Anthony Dominguez, Jan Kleissl, Janet Shields, and Byron Washom, Intra-hour forecasting with a total sky imager at the UC San Diego solar energy testbed, Solar Energy 85 (2011), no. 11, 2881–2893. [6] Yinghao Chu, Hugo T.C. Pedro, and Carlos F.M. Coimbra, Hybrid intra-hour DNI forecasts with sky image processing enhanced by stochastic learning, Solar Energy (2013). [7] Soumyabrata Dev, Yee Hui Lee, and Stefan Winkler, Color-Based Segmentation of Sky/Cloud Images from Ground-Based Cameras, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 10 (2017), no. 1, 231–242. [8] H Escrig, F J Batlles, J Alonso, F M Baena, J L Bosch, I B Salbidegoitia, and J I Burgaleta, Cloud detection, classification and motion estimation using geostationary satellite imagery for cloud cover forecast, 2013. [9] Hao Huang, Shinjae Yoo, Dantong Yu, Dong Huang and Hong Qin, Cloud Motion Detection for Short Term Solar Power Prediction, ICML 2011 Workshop on Machine Learning for Global Challenges, 2011. [10] Qingyong Li, Weitao Lu, and Jun Yang, A hybrid thresholding algorithm for cloud detection on groundbased color images, Journal of Atmospheric and Oceanic Technology 28 (2011), no. 10, 1286–1296. [11] Qingyong Li, Weitao Lu, Jun Yang, and James Z. Wang, Thin cloud detection of all-sky images using Markov random fields, IEEE Geoscience and Remote Sensing Letters 9(2012), no. 3, 417–421. [12] Charles N. Long, J. M. Sabburg, J. Calbó, and D. Pagès, Retrieving cloud characteristics from groundbased daytime color all-sky images, Journal of Atmospheric and Oceanic Technology (2006). [13] Bruce D. Lucas and Takeo Kanade, An Iterative Image Registration Technique with an Application to Stereo Vision., Robotics (1981). [14] Lydie Magnone, Fabrizio Sossan, Enrica Scolari, and Mario Paolone, Cloud Motion Identification Algorithms Based on All-Sky Images to Support Solar Irradiance Forecast. 45