Monitorización y Predicción de la Radiación Solar Mediante Visión del Paso de Nubes

Monitorizaci´on y Predicci´on de la Radiaci´on Solar Mediante Visi´on del Paso de Nubes Alberto Mo˜nino Hormigo Alb[email protected] Francisco Rodr´ıguez Rubio [email protected] Manuel Ruiz Arahal [email protected] Resumen En este art´ıculo se expone un m´etodo para conocer qu´e fracci´on de un campo solar se encontrar´a sombreada por el efecto de una nube interpuesta en la linea que une el campo solar con el sol mediante el an´alisis de una secuencia de im´agenes tomadas antes de que suceda dicho acontecimiento. Para obtener la posici´on futura de una nube es necesario calcular su velocidad. Para ello se utiliza el m´etodo de Lucas-Kanade, que analiza el flujo ´optico para obtener la velocidad de algunos de los puntos que se encuentran dentro de la matriz de inercia equivalente de una nube. Para obtener la velocidad aproximada de la nube se realiza una media de velocidades obtenidas de dichos puntos. Para realizar c´alculos en la imagen, como las im´agenes tomadas son hemisf´ericas, ha de realizarse un proceso para la eliminaci´on de la distorsi´on, as´ı como una serie de tratamientos de imagen.Por ´ultimo, si se supone que las nubes se encuentran todas a una misma altura y que no tienen fondo se puede calcular la proyecci´on de la nube predicha sobre el campo solar con una proyecci´on paralela y comprobar si los puntos del campo solar pertenecen a la elipse proyectada predicha, obteniendo as´ı qu´e fracci´on del campo se encontrar´a oculta bajo una nube. Palabras clave: Predicci´on, Paso de nubes, Despecificaci´on, M´etodo de Lucas-Kanade, Momentos Invariantes. 1. Introducci´on Con el gran auge de las energ´ıas renovables en la actualidad, surgen nuevas formas para la optimizaci´on de la producci´on de dicha energ´ıa. En concreto, en la producci´on de energ´ıa a partir de la transformaci´on de la radiaci´on solar se est´a investigando actualmente en nuevas formas de control de los campos solares basadas en la predicci´on de la cobertura de nubes. Espa˜na es actualmente uno de los pa´ıses pioneros en la investigaci´on e implantaci´on de este tipo de instalaciones para la obtenci´on de energ´ıas limpias. Espa˜na posee una localizaci´on muy adecuada para la instalaci´on de este tipo de tecnolog´ıa, ya que la irradiancia media recibida a lo largo del a˜no en nuestro pa´ıs es muy elevada. Para el control de un campo solar, tradicionalmente se han usado modelos de alta resoluci´on no lineales y con una gran din´amica. Para obtener predicciones en base a estos modelos se utilizaban m´etodos num´ericos as´ı como por im´agenes obtenidas a partir de sat´elites. Estas predicciones de la radiaci´on solar tienen un horizonte de predicci´on bastante grande, de alrededor de 30 minutos, y son utilizados para grandes superficies, ya que, aunque un sat´elite sea un excelente m´etodo de obtener informaci´on en estos casos, no tiene la suficiente precisi´on para regiones peque˜nas, lo que hace que la localizaci´on exacta de una nube resulte dif´ıcil. Estos modelos est´an muy expuestos a perturbaciones producidas en cortos espacios de tiempo, situaci´on muy habitual en los sistemas meteorol´ogicos, por lo que hace realmente dif´ıcil obtener una predicci´on precisa a corto plazo. [2]. Debido a estos motivos, surge la necesidad de buscar nuevos sistemas de adquisici´on de datos con mayor precisi´on, que obtengan informaci´on local del sistema. En este contexto surge la utilizaci´on de c´amaras que captan im´agenes de las nubes como m´etodo de adquisici´on de datos locales para predicciones a corto plazo, as´ı como de redes de sensores de radiaci´on. Las redes de sensores precisan de grandes superficies para instalar todos los sensores a una distancia suficientemente alejada del campo para obtener predicciones de la cobertura de nubes y no datos en tiempo real y con una separaci´on adecuada entre ellos para que no fuesen datos redundantes. Esto supondr´ıa que habr´ıa que poseer grandes superficies cercanas al campo para su instalaci´on, lo que en muchos casos hace inviable su utilizaci´on. En contraposici´on, con pocas c´amaras instaladas alrededor del campo solar es posible obtener una caracterizaci´on precisa de la cobertura de nubes actual. Adem´as, mediante otras t´ecnicas basadas en el flujo ´optico es posible obtener predicciones de la distribuci´on tridimensional de las nubes a cor- to plazo analizando secuencias de im´agenes. Para obtener estimaciones tridimensionales es impres- Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 91 cindible la utilizaci´on de al menos dos c´amaras, pero sigue siendo un n´umero mucho menor que el que es necesario para construir una red de sensores, adem´as de que con una red de sensores no es posible obtener una estimaci´on tridimensional. Por este motivo, el uso de este tipo de dispositivos para la predicci´on de la radiaci´on solar a corto plazo es muy prometedor. Para la elaboraci´on de este trabajo no se ha realizado el c´alculo de la altura, s´olo se ha obtenido la posici´on de la nube en un plano paralelo al suelo. Se supondr´a a partir de ahora una altura conocida para la realizaci´on de las tareas. 2. Adquisici´on de im´agenes La c´amara que se ha utilizado para la toma de im´agenes es la Q-24M de Mobotix. Se trata de una c´amara hemisf´erica con una resoluci´on de hasta 2048x1536, de bajo coste en comparaci´on con las c´amaras que se utilizan normalmente para tomar im´agenes del cielo. Las im´agenes adquiridas son RGB. Para controlar la c´amara, ´esta dispone de un software instalado en su memoria, al cual se accede desde el computador a trav´es del navegador. Para el almacenamiento de las im´agenes se programa en la c´amara un evento peri´odico que toma una imagen cada 5 segundos y la transfiere al ordenador mediante un protocolo SMB(Server Message Block). 3. Tratamiento de im´agenes Se han realizado una serie de tareas para poder mejorar el posterior an´alisis y extracci´on de informaci´on de la imagen. Las tareas que se describen en esta secci´on han de aplicarse a toda la secuencia de im´agenes que se necesitan para el c´alculo de velocidades, realizado en una secci´on posterior. 3.1. Despecificaci´on En primer lugar se ha realizado una transformaci´on llamada despecificaci´on, cuyo objetivo es solucionar la distorsi´on introducida en la captaci´on de la imagen por ser la c´amara hemisf´erica. La imagen obtenida por la c´amara tiene una funci´on de mapeo ortogr´afico. La funci´on de transformaci´on que se necesita para remapear la imagen a una proyecci´on rectil´ınea es: rop =fop ∗sin(θ) (1) siendo rop la distancia de un punto de la imagen original al centro de ´esta, fop la distancia focal de la c´amara hemisf´erica y θel ´angulo en radianes que forma el punto con el eje ´optico. El valor de la distancia focal de la c´amara no es un dato, por tanto, se elige un valor que haga que la despecificaci´on sea adecuada por prueba y error. Se recorta la parte central de la imagen, puesto que las nubes que se muestran fuera de esa zona est´an a una distancia muy elevada, por lo que no ser´a hasta pasado bastante tiempo cuando se vean, y ya se habr´an tomado otras im´agenes cuando eso ocurra, ademas de que esa zona de la imagen se encuentra bastante deformada. Esto reduce la carga computacional ya que son muchos menos p´ıxeles los que han de ser transformados. Ahora se realiza la transformaci´on para cada pixel, excepto para el centro de la imagen, ya que tiene los mismos valores que el pixel de la imagen original. En primer lugar se calcula el radio del pixel en la imagen original, es decir, la distancia al centro de la imagen y el ´angulo seg´un las siguientes f´ormulas: r=q(xpixel −xcentro)2+ (ypixel −ycentro)2 (2) θ= arctan(xcentro −xpixel ycentro −ypixel ) (3) El ´angulo θse mantiene constante en la transformaci´on, por lo que tan solo falta por calcular el radio que tendr´a el pixel en la imagen despecificada. En la imagen despecificada, el radio de la imagen se calcular´a de la siguiente forma: r=f∗tan(θ) (4) La distancia focal que se utiliza en la ecuaci´on (4) no es dato, por lo que, como sucede con la distancia focal de la imagen hemisf´erica, se elige un valor por prueba y error que haga que la despecificaci´on sea adecuada. Para calcular el radio de la imagen despecificada basta con despejar θde la ecuaci´on (1) y sustituirlo en la ecuaci´on (4) quedando la ecuaci´on: r=f∗tan(arcsin( rop fop )) (5) Una vez que se tiene el radio y el ´angulo del pixel, basta con ir al pixel original, que est´a totalmente definido, y copiar su intensidad al pixel de la imagen despecificada. Este tratamiento se realiza s´olo para dos de los planos de color de la imagen, que ser´an el verde y el rojo, tomando 0.6 de rojo y 0.4 de verde. El azul no se utiliza pues es el color del fondo, es decir, del cielo. A continuaci´on se compone la nueva imagen con estos dos planos de colores. Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 92 3.2. Transformaci´on polin´omica para corregir la distorsi´on radial en la lente Tras la despecificaci´on se obtiene una imagen que sigue teniendo una distorsi´on de tipo coj´ın (pincushion distortion) producida por la lente. Se produce porque los puntos de la imagen se desplazan respecto a lo que indica la ecuaci´on de proyecci´on de la c´amara estenopeica, sobre todo la componente radial. Esto quiere decir que en vez de proyectarse sobre una linea recta, la lente distorsiona el ´angulo con el que llegan los fotones, produciendose el desplazamiento antes comentado. Para resolver este problema es necesario aplicar a la imagen una transformaci´on polin´omica para corregir la distorsi´on radial en la lente. La transformaci´on polin´omica se aplica a todos los p´ıxeles de la imagen y sigue la siguiente ecuaci´on: (xf, yf)=(xi, yi)∗(1 + k1∗ri) (6) siendo (xf, yf) las coordenadas del punto tras la transformaci´on, (xi, yi) las del punto inicial, y ri el radio del pixel (xi, yi). La variable k1se ha obtenido por prueba y error, sabiendo que tiene que ser mayor que cero, y que la imagen debe reducir su tama˜no. El algoritmo recorre la nueva imagen copiando los p´ıxeles equivalentes de la imagen distorsionada. Se podr´ıan incluir mas t´erminos en el polinomio que es funci´on del radio, pero no es necesario ya que no se busca una transformaci´on con gran precisi´on. El problema de este m´etodo es que el sistema de ecuaciones (5) es un sistema incompatible. if=ii∗(1 + k1∗ri) jf=ji∗(1 + k1∗ri) ri=pi2 i+j2 i (7) Para encontrar una soluci´on a este problema se crea una tabla de equivalencias entre radios partiendo de la ecuaci´on (6). Ser´an dos vectores con el mismo n´umero de componentes en los que el valor de una posici´on corresponde con la transformaci´on del pixel que ocupa la misma posici´on en el otro vector. rf=qi2 f+j2 f =pi2 i∗(1 + k1∗ri)2+j2 i∗(1 + k1∗ri)2 = (1 + k1∗ri)∗pi2 i+j2 i=ri∗(1 + k1∗ri) (8) El procedimiento consiste en recorrer la nueva imagen a completar, calculando el radio que tiene cada pixel, buscando en el vector de radios transformados el radio m´as cercano al obtenido, y se obtiene el equivalente con distorsi´on sabiendo que ocupar´a la misma posici´on en el vector de radios sin distorsi´on y se calcula el ´angulo θ, que es igual en ambas im´agenes, utilizando la ecuaci´on (9). θ= arctan( if jf ) (9) Teniendo el radio y el ´angulo de los p´ıxeles de ambas im´agenes, se conocen las coordenadas del pixel de la imagen distorsionada y de la nueva imagen, entonces se copia la intensidad del pixel de la imagen distorsionada a la nueva imagen, ya sin distor- si´on. Por ´ultimo, a la imagen obtenida tras la despecificaci´on y la transformaci´on polin´omica hay que realizarle un suavizado. Este proceso es necesario para solucionar el problema debido a la aproximaci´on realizada durante la transformaci´on polin´omica, tras la cual, algunos p´ıxeles de la imagen no fueron transformados, quedando sin intensidad. 4. Caracterizaci´on de nubes Para obtener las propiedades de las nubes se ha de binarizar la imagen obtenida tras el tratamiento. La binarizaci´on se realiza con un umbral variable. La primera propiedad necesaria es el ´area de las regiones. ´ Esta se utilizar´a para eliminar aquellas regiones que sean peque˜nas o grandes. Ser´an peque˜nas aquellos puntos que no fueron eliminados con el suavizado, as´ı como aquellas nubes que son muy peque˜nas y no son necesaria tenerlas en cuenta. Las grandes el destello producido por el sol, que en caso de que ninguna nube se encuentre muy pr´oxima de forma que disminuya la intensidad de los rayos, tendr´a un tama˜no mucho mayor que las nubes. Para eliminar dichas regiones, basta con vaciar las posiciones que ocupan las regiones en el vector de estructuras. Las otras propiedades ser´an el centro de gravedad, los ejes de la elipse de inercia equivalente y la orientaci´on de la elipse. Las elipses se utilizar´an para aproximar la forma de las nubes. Con estas propiedades se tiene caracterizada la posici´on y el tama˜no de una nube en la imagen. 5. C´alculo de velocidades La ´ultima caracter´ıstica a extraer de las im´agenes es la velocidad de una nube a partir de una secuencia de im´agenes tomadas en un ∆t, con una separaci´on de tiempo entre im´agenes constante.El objetivo de este c´alculo es poder estimar la posici´on de la nube pasado un tiempo ∆t0. Esto es sencillo, pues conociendo el ∆t0, la velocidad de un punto determinado y su posici´on. Para obtener la velocidad se va a utilizar el flujo ´optico, que Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 93 es el patr´on de velocidades aparentes que tienen lo objetos con respecto de la c´amara. Se parte de una secuencia de im´agenes, las cuales tienen p´ıxeles con un nivel de gris I(x, y). Como son varias las im´agenes y se encuentran ordenadas en el tiempo se a˜nade una nueva coordenada t, que es la posici´on de la imagen dentro de la secuencia, teniendo as´ı I(x, y, t). Si se supone que la luminosidad no var´ıa, para estimar el flujo ´optico se puede considerar que el nivel de gris de un objeto no va a variar entre dos im´agenes consecutivas de la secuencia al producirse el movimiento, ya que la velocidad de las nubes es peque˜na en comparaci´on con el tiempo entre una imagen y otra. ∂I ∂xVx+∂I ∂y Vy+∂I ∂t = 0 (10) En la ecuaci´on (10) se tiene una ecuaci´on con dos inc´ognitas, por lo que existen infinitas soluciones para resolverse. Para hallar una soluci´on a este problema se aplica el M´etodo de Lucas-Kanade. Se acepta la hip´otesis de que p´ıxeles cercanos dentro de la imagen poseen la misma velocidad, por lo que aplicando la misma ecuaci´on a distintos p´ıxeles con las mismas inc´ognitas se obtiene un sistema sobredeterminado, que ya si es posible resolverlo , por ejemplo, mediante el m´etodo de los m´ınimos cuadrados, obteniendo una estimaci´on de la velocidad. Los t´erminos que tienen derivadas parciales se pueden reescribir como gradientes, obteniendo el siguiente sistema de ecuaciones para los puntos p, q, r... pertenecientes a la vecindad considerada: Ix(p)∗Vx+Iy(p)∗Vy=−It(p) Ix(q)∗Vx+Iy(q)∗Vy=−It(q) Ix(r)∗Vx+Iy(r)∗Vy=−It(r) . . . (11) Este sistema de ecuaciones puede reescribirse de forma matricial de la forma A∗v=bcon A=     Ix(p)Iy(p) Ix(q)Iy(q) Ix(r)Iy(r) . . .. . .      ,b=−      It(p) It(q) It(r) . . .      yv=Vx Vy(12) Para obtener la velocidad de la nube se va a obtener una media de velocidades. Se calcular´an las velocidades aproximadas de 20 vecindades situadas sobre los ejes de la elipse que representa la nube. Este c´alculo se realiza para todas las nubes detectadas en la tarea anterior. El resultado del M´etodo de Lucas-Kanade ser´a la velocidad media de la nube aplicada al centro de gravedad de la nube. Teniendo una secuencia de m´as de dos im´agenes se pueden calcular m´as velocidades de esa nube a lo largo del tiempo siempre y cuando la nube no desaparezca. Con todas las velocidades obtenidas de la secuencia se vuelve a realizar una media entre todas ellas y una vez que se tiene esa velocidad es posible calcular la posici´on de la elipse tras un tiempo ∆t0aplicando la ecuaci´on (13) que se muestra a continuaci´on. (xf, yf)=(xi, yi) + vm∗∆t0(13) 6. Obtenci´on de coordenadas 3D Se pretende obtener los par´ametros intr´ınsecos de la c´amara para la obtenci´on de las coordenadas 3D de las nubes aplicando las ecuaciones de proyecci´on de una c´amara estenopeica. En primer lugar, las posiciones de las nubes que se han obtenido est´an referenciadas a la esquina superior izquierda de la imagen. Para trabajar con distancias reales es necesario trasladar este sistema de coordenadas a uno con origen en el centro de la imagen, es decir, a un sistema de coordenadas con origen en lo que ser´ıa el estenopo en una c´amara estenopeica, y a continuaci´on, transformar las distancias obtenidas en p´ıxeles a metros, ya que el campo solar est´a medido en metros. Para trasladar el sistema de referencia a la posici´on de la c´amara hay que restar a las coordenadas de todos los centros de gravedad de las elipses la mitad del n´umero de filas o columnas de la imagen.Adem´as, como las im´agenes fueron tomadas desde el suelo, con un sistema de coordenadas en Y que apuntaba al sur y un eje X que apuntaba al oeste, para realizar las proyecciones sobre el suelo, deben estar en un sistema de coordenadas coherente. En el sistema de referencia del mundo el eje Y apunta hacia el norte y el eje X hacia el este. Es necesario realizar una rotaci´on de 180 grados de todas las coordenadas obtenidas. Para ello se multiplican las coordenadas obtenidas por la matriz de rotaci´on de la ecuaci´on (14). R=cos(180) sin(180) sin(180) cos(180) (14) Para hacer la transformaci´on de p´ıxeles a metros hay que hacer uso de la geometr´ıa proyectiva. Las ecuaciones de proyecci´on de un punto en el espacio 3D a un punto en la imagen tienen la siguiente forma: x=−f∗X Zy=−f∗Y Z(15) En lugar de la distancia focal en metros, se van a transformar las ecuaciones (15) para que no dependan ni de la posici´on del punto en la imagen, Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 94 que est´a medida en metros, ni de la distancia focal, tambi´en medida en metros. Aplicando las ecuaciones (16) y (17) es posible obtener dichas ecuaciones pero dependiendo de par´ametros que est´en medidos en p´ıxeles. Se sabe que la distancia focal en p´ıxeles viene dada por la ecuaci´on (16) donde mxes un par´ametro de la c´amara que se interpreta como un factor de conversi´on de p´ıxeles a metros. Si se aplica este factor de conversi´on a una distancia real del mundo 3D, por ejemplo, a la del eje X se obtiene la ecuaci´on (17) donde ∆jes una distancia en la imagen medida en p´ıxeles, y xes una distancia medida en metros. αx=f∗mx(16) x=∆j mx (17) Si se introduce las ecuaciones (17) y (16) en la ecuaci´on (15) se puede despejar αxyαy, pues todos los dem´as t´erminos son conocidos en las im´agenes que se est´an analizando, quedando las ecuaciones que se muestran en la ecuaci´on (18). Ya no se va a obtener fde las im´agenes, sino que se va a obtener αxyαy. αx=−∆j∗Z Xαy=−∆i∗Z Y(18) Para obtener estos par´ametros se ha utilizado una im´agen tomada de una cuadricula en la que se puede medir la distancia entre dos puntos y la distancia a la que se encuentra la cuadr´ıcula de la c´amara. Las ecuaciones 18, para el caso que se est´a tratando en el que Z es constante, se pueden utilizar indistintamente para puntos como para distancias, pues al aplicar incrementos, se obtienen incrementos. Con esta imagen, de las que se conocen propiedades de objetos que aparecen en ella, se eliminan las inc´ognitas Z,XeY, quedando por tanto s´olo las inc´ognitas αxyαy. La imagen utilizada para calcular estos par´ametros se muestran en la figura 1 De la imagen se obtiene un valor de αxyαyigual para ambos de 163. Finalmente, reordenando las ecuaciones (18), y una vez ya se han obtenido los par´ametros αxyαyse puede obtener las coordenadas de un punto en el espacio 3D con las ecuaciones (19) en funci´on de las coordenadas del punto en la imagen medido en p´ıxeles con respecto al centro de la imagen. X=−j∗Z αx Y=−i∗Z αy(19) Se puede ver que estas ecuaciones dependen de la distancia del objeto a localizar en el mundo real en Figura 1: Imagen cuadriculada para la obtenci´on de los par´ametros intr´ınsecos de la c´amara. el eje Z. Por este motivo, las transformaciones inversas desde el mundo 2D al mundo 3D no siempre son posibles. Es necesario obtener este par´ametro mediante otro m´etodo. Una posible soluci´on que no est´a contemplada en este trabajo es usar m´as c´amaras para realizar una triangulaci´on del pun- to 3D y as´ı obtener la coordenada Z del punto. Obteniendo esta coordenada basta con aplicar las ecuaciones (19) para obtener las otras dos coordenadas del punto y estar as´ı totalmente definido. Los semiejes de las elipses obtenidos de las im´agenes est´an medidos en p´ıxeles. Para trabajar con ellos es necesario convertirlos a metros. Como en este caso αxyαytienen el mismo valor, se puede aplicar directamente las ecuaciones 19 para incrementos de x o de y indistintamente. Por tanto, para obtener los valores en metros basta con sustituir joipor un ∆jo ∆i, y se obtendr´a ∆Xo ∆Yrespectivamente, medido en metros. Para la orientaci´on no es necesario realizar ning´un tipo de conversi´on. 7. C´alculo de proyecciones En esta ´ultima secci´on se pretende obtener qu´e fracci´on del campo solar se ver´a afectada por la nube tras un ∆tpara que el algoritmo de control pueda anticipar estas perturbaciones. Los datos de entrada son el centro de gravedad de la elipse predicha, sus ejes y su orientaci´on. A partir de ahora se supondr´a que las nubes se encuentran en un plano paralelo al suelo, ya que no es posible con una ´unica c´amara conocer la profundidad de una nube, imposibilitando la estimaci´on del volumen de la nube. Como la altura no se ha podido calcular, se supondr´a una altura Z=500 para todas las nubes quedando como un par´ametro para que en el futuro se pueda incorporar como variable de entrada. En primer lugar, para conocer la proyecci´on de la Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 95 nube es necesario conocer la fuente de luz que provoca dicha proyecci´on, que en este caso es el sol. Para obtener la posici´on del sol mediante las ecuaciones solares se utiliza un algoritmo presentado en [1] e implementado en MATLAB por Vincent Roy. Consiste en una funci´on con argumento de entrada la posici´on de la c´amara y la hora en la que se tom´o la imagen y con argumento de salida el zenit y acimut con respecto a la posici´on de la c´amara, que definen completamente la posici´on del sol. Una vez se tiene la posici´on del sol, se procede al c´alculo de las proyecciones. Se puede suponer que el sol se encuentra a una distancia infinita, por tanto, los rayos procedentes de ´el vienen desde el infinito, simplificando esto mucho los c´alculos, ya que las lineas procedentes de un punto en el infinito son todas paralelas. Se realiza entonces una proyecci´on paralela. Como se ha supuesto que las nubes se encuentran en un plano paralelo al suelo, s´olo es necesario proyectar el centro de gravedad de la elipse, ya que los ejes de la elipse y su orientaci´on no cambian al realizar la proyecci´on. Para calcular la posici´on del centro gravedad proyectado basta con calcular la ecuaci´on de la recta que pasa por el centro de gravedad a una altura Z y que tiene la misma direcci´on que el vector de direcci´on de los rayos solares, que est´a totalmente definido por las ecuaciones solares. Siendo θel zenit, que se mide desde el suelo y ψel acimut, medido desde el norte con sentido positivo hacia el este se tiene que: i= sin(π−ψ) (20) j=−cos(π−ψ) (21) k= tan(θ) (22) con ien direcci´on este y jen direcci´on norte. A partir de la ecuaci´on (20) el vector de direcci´on de los rayos de sol tendr´a la misma direcci´on que (i,j,k), pero con sentido opuesto. A partir de las coordenadas del centro de gravedad de la elipse (xe, ye, ze) y del vector de direcci´on u se obtiene la ecuaci´on de la recta buscada. x−xe −sin(π−ψ)=y−ye cos(π−ψ)=z−ze −tan(θ)(23) Si se conoce que la cota del punto buscado es z= 0, se puede despejar xpeyp, coordenadas del punto proyectado, de un sistema de ecuaciones de dos ecuaciones y dos inc´ognitas, que posee una ´unica soluci´on. xp=xe−ze∗sin(π−ψ) tan(θ) yp=ye+ze∗cos(π−ψ) tan(θ) (24) El algoritmo pensado para conocer la fracci´on del campo que se encontrar´a perturbado por el efec- to de la nube consiste en conocer qu´e puntos del campo solar se encuentran dentro de la elipse que simula su efecto proyectada. Un punto del campo cuya suma de distancias a los focos de la elipse proyectada sea menor al tama˜no del eje mayor de dicha elipse se encontrar´a sujeta a la perturbaci´on de una nube, es decir, que estar´a sombreada, siendo ael semieje mayor. si P F1+P F2<2∗a=⇒P∈elipse (25) La posici´on de los focos vendr´a dada por las ecuaciones (26) donde αdefine la orientaci´on de la nube, (xF1, yF1)y(xF2, yF2) la posici´on de los focos 1 y 2 respectivamente y (xp, yp) la posici´on de la proyecci´on del centro de gravedad y c=√a2−b2. xF1=xp+c∗cos(α) yF1=yp+c∗sin(α) xF2=xp−c∗cos(α) yF2=yp−c∗sin(α) (26) Por ´ultimo, se comprueba la ecuaci´on (25) para todos los puntos del campo solar, completando una matriz de ceros del mismo tama˜no que el campo con un uno en caso de que se cumpla. Como se ha supuesto que el centro de coordenadas se encuentra en la c´amara, pues es en el sistema de referencia en el que se calcul´o el centro de la elipse, para hacer la comprobaci´on es necesario a˜nadir la posici´on del campo solar con respecto al sistema de referencia, es decir, con respecto a la c´amara. Suponiendo que el campo solar sea rectangular y que sea paralelo a los ejes de coordenadas, si Lxes la longitud del campo en xyLyes la longitud del campo en y, y la posici´on de la esquina superior izquierda del campo se encuentra en (xcampo, ycampo), la condici´on se puede formular de la siguiente forma: j∈[1, Lx] i∈[1, Ly] H(i, j) = 1 si p(xcampo +j−xF1)2+ (ycampo −i−yF1)2+ p(xcampo +j−xF2)2+ (ycampo −i−yF2)2<2∗a (27) donde H(i, j) es la matriz compuesta por unos y ceros que indica que partes del campo solar se predice que estar´an cubiertas tras un tiempo ∆t0in- troducido por el operario, es decir, que es la entrada del algoritmo de control que se buscaba. 8. Resultados Se muestran los resultados obtenidos con los algoritmos rese˜nados anteriormente. Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 96 En primer lugar, en la figura 2 se muestra una imagen original, antes de realizar ning´un acci´on sobre ella. A continuaci´on, en la figura 3 se muestra la imagen obtenida tras realizar el tratamiento de la imagen. Figura 2: Imagen original. Figura 3: Imagen tratada. En la figura 4 se pueden observar los resultados obtenidos de una secuencia de 7 im´agenes, con un intervalo entre ellas de 5 segundos. La X verde marca la posici´on predicha de la imagen tras un ∆tde 5 minutos y 40 segundos. Las elipses que aparecen tienen un color m´as negro cuanto m´as antigua sea la imagen. Se puede observar como la posici´on obtenida es coherente con el movimiento que tiene la nube. La precisi´on de la posici´on estimada se puede comprobar en la figura 5. Se puede observar, comparando con la posici´on predicha tras el an´alisis, que el m´etodo de predicci´on utilizado tiene una preci- si´on considerable, aunque se a˜nade que con horizontes de predicci´on elevados, el m´etodo est´a muy expuesto a perturbaciones del tipo rachas de vien- to. La proyecciones en el campo se pueden observar en la figura 6. En el caso que se muestra, ninguna nube proyecta su sombra sobre el campo, por Figura 4: Posici´on estimada tras el an´alisis de una secuencia. Figura 5: Imagen para comprobar los resultados. Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 97 lo que la matriz H resultante ser´a una matriz de ceros. Figura 6: Posici´on estimada tras el an´alisis de una secuencia. 9. Conclusi´on Es posible obtenerse proyecciones precisas de las nubes para condiciones en la que no existen perturbaciones como es el viento. Cuando se habla del viento como perturbaci´on se refiere a que exista una variaci´on considerable en su velocidad o direcci´on. Al ser el algoritmo dise˜nado para cualquiera que sea el n´umero de nubes que haya en el momento, da cierta versatilidad, ya que no hay que centrarse en una ´unica nube, sino que pueden ser estudiadas todas a la vez, sin olvidar que las nubes que se mueven en una misma imagen tienen una relaci´on, ya que aproximadamente todas las nubes se mueven con la misma direcci´on. La velocidad tambi´en debe ser parecida en todas las nubes, pues si no existen perturbaciones, todas est´an sometidas a las mismas corrientes de viento. Esto permite realizar ciertas comprobaciones de la veracidad de los resultados obtenidos. Las condiciones meteorol´ogicas son muy variables a lo largo del tiempo, por lo que las predicciones que se realizan est´an muy expuestas a las perturbaciones antes comentadas. Por este motivo, se recomienda que los horizontes de predicci´on utilizados no sean excesivamente largos. Tambi´en se recomienda que las c´amaras no se sit´uen muy alejadas del campo solar, pues esto conllevar´ıa a realizar predicciones con grandes horizontes de predicci´on. Se proponen aqu´ı, como soluci´on a tomar para mejorar el problema del sol, tomar tres im´agenes en vez de una, como se explica en Solar Energy Forecasting and Resource Assessement [2], y realizar una combinaci´on de las tres formando lo que se llama una imagen HDR. En este libro se propone tomar im´agenes con tiempos de exposici´on diferentes, que son de 5, 20, y 80 milisegundos. Se necesita, por tanto, una c´amara que permita realizar un ajuste manual del tiempo de exposici´on. Es necesaria la instalaci´on de m´as de una c´amara para obtener la altura de las nubes. Esta caracter´ıstica es totalmente necesaria para la estimaci´on del porcentaje de campo que se encuentra tapado por la nube. Para hallar la altura de una nube se podr´ıa realizar una triangulaci´on de puntos singulares, como podr´ıan ser los extremos de las nubes, pues si se conoce la posici´on del punto en las distintas im´agenes, y se conoce la posici´on relativa de las c´amaras la obtenci´on es trivial resolviendo las distintas ecuaciones obtenidas. En caso de que se obtuvieran varias alturas debido a la profundidad que las nubes poseen, se podr´ıa utilizar la menor, pues es la que afecta en mayor medida a las proyecciones. La instalaci´on de m´as c´amaras tambi´en puede ayudar a comprobar resultados, en caso de que no se encuentren muy alejadas y sean capaces de observar la misma nube. A pesar de ello, Este c´alculo forma parte de un trabajo diferente, que no afecta a los algoritmos expuestos en este art´ıculo. Agradecimientos Los autores desean expresar su agradecimiento al ministerio de Ciencia e Innovaci´on por la financiaci´on de este trabajo con el proyecto de referencia DPI2013-44135-R. Referencias [1] Ibrahim Reda and Afshin Andreas, (2008) Solar Position Algorithm for Solar Radiation Applications, National Renewable Energy Laboratory (NREL). [2] Jan Kleissl, (2013)Solar Energy Forecasting and Resource Assessement, pp 196- 228,Academic Press. [3] Javier Gonzalez-Jimenez, (1999) Visi´on por Computador, Thomson-Paraninfo. [4] Manuel Ruiz Arahal, Apuntes de Sistemas de Percepci´on. [5] Mobotix, Camera Manual Q24M. [6] Wikipedia contributors(25 May 2015) Fisheye lens, Wikipedia, The Free Encyclopedia, 663995745. Actas de las XXXVI Jornadas de Automática, 2 - 4 de septiembre de 2015. Bilbao ISBN 978-84-15914-12-9 © 2015 Comité Español de Automática de la IFAC (CEA-IFAC) 98