FR3047829A1 - METHOD FOR ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE - Google Patents
METHOD FOR ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE Download PDFInfo
- Publication number
- FR3047829A1 FR3047829A1 FR1651122A FR1651122A FR3047829A1 FR 3047829 A1 FR3047829 A1 FR 3047829A1 FR 1651122 A FR1651122 A FR 1651122A FR 1651122 A FR1651122 A FR 1651122A FR 3047829 A1 FR3047829 A1 FR 3047829A1
- Authority
- FR
- France
- Prior art keywords
- image
- inn
- digital
- calculating
- solar disk
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 84
- 238000003384 imaging method Methods 0.000 claims abstract description 24
- 230000003287 optical effect Effects 0.000 claims abstract description 23
- 239000013598 vector Substances 0.000 claims abstract description 19
- 238000005457 optimization Methods 0.000 claims abstract description 8
- 238000009826 distribution Methods 0.000 claims description 64
- 239000011159 matrix material Substances 0.000 claims description 43
- 238000004364 calculation method Methods 0.000 claims description 38
- 238000006073 displacement reaction Methods 0.000 claims description 19
- 238000004458 analytical method Methods 0.000 claims description 14
- 230000007717 exclusion Effects 0.000 claims description 14
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 11
- 125000001475 halogen functional group Chemical group 0.000 claims description 8
- 230000004907 flux Effects 0.000 claims description 6
- 230000010365 information processing Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 description 18
- 230000006870 function Effects 0.000 description 17
- 238000004519 manufacturing process Methods 0.000 description 17
- 238000001514 detection method Methods 0.000 description 15
- 230000011218 segmentation Effects 0.000 description 14
- 238000011282 treatment Methods 0.000 description 14
- 230000005855 radiation Effects 0.000 description 8
- 238000013459 approach Methods 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 6
- 238000012937 correction Methods 0.000 description 6
- 230000002123 temporal effect Effects 0.000 description 6
- 238000012360 testing method Methods 0.000 description 6
- 230000005611 electricity Effects 0.000 description 5
- 230000015654 memory Effects 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 3
- 230000005484 gravity Effects 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000007781 pre-processing Methods 0.000 description 3
- 238000011084 recovery Methods 0.000 description 3
- 230000004075 alteration Effects 0.000 description 2
- 230000003466 anti-cipated effect Effects 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 230000001960 triggered effect Effects 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 239000011248 coating agent Substances 0.000 description 1
- 238000000576 coating method Methods 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000002513 implantation Methods 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000001681 protective effect Effects 0.000 description 1
- 230000035484 reaction time Effects 0.000 description 1
- 230000011514 reflex Effects 0.000 description 1
- 230000000284 resting effect Effects 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/90—Determination of colour characteristics
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
- G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
- G06T7/75—Determining position or orientation of objects or cameras using feature-based methods involving models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10016—Video; Image sequence
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10024—Color image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30192—Weather; Meteorology
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30244—Camera pose
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
- Length Measuring Devices By Optical Means (AREA)
Abstract
Un procédé d'estimation de la position du disque solaire dans une image numérique de ciel acquise par un dispositif d'acquisition comprenant une optique et un capteur numérique disposé dans le plan focal de l'optique, l'estimation de la position du disque étant réalisée en fonction d'un modèle astronomique de l'azimut et de l'élévation du soleil dans le ciel ayant pour paramètre une position prédéterminée (θr1, θr2,Azr) du dispositif d'imagerie. Le procédé comprend une étape de détermination de la position prédéterminée (θr1, θr2,Azr) du dispositif d'imagerie, comportant : - l'acquisition d'un flux d'images numériques de calibration horodatées {INcn} comprenant chacune le disque solaire ; - pour chaque image INcn dudit flux {INcn}, le calcul d'une position Pcn du centre du disque solaire dans l'image INcn en analysant l'image INcn ; - le calcul de la position du dispositif d'imagerie (θr1, θr2, Azr) comme solution d'un problème d'optimisation selon la relation : où Ω est un ensemble de recherche, et || || est une fonction de coût quantifiant un écart entre deux vecteurs de R3.A method of estimating the position of the solar disk in a digital sky image acquired by an acquisition device comprising an optical and a digital sensor disposed in the focal plane of the optics, the estimation of the position of the disk being performed according to an astronomical model of the azimuth and the elevation of the sun in the sky having for parameter a predetermined position (θr1, θr2, Azr) of the imaging device. The method comprises a step of determining the predetermined position (θr1, θr2, Azr) of the imaging device, comprising: - acquiring a stream of calibrated digital calibration images {INcn} each comprising the solar disk; for each image INcn of said stream {INcn}, calculating a position Pcn of the center of the solar disk in the image INcn by analyzing the image INcn; calculating the position of the imaging device (θr1, θr2, Azr) as a solution of an optimization problem according to the equation: where Ω is a search set, and || || is a cost function quantifying a difference between two vectors of R3.
Description
PROCEDE D’ESTIMATION DE LA POSITION DU DISQUE SOLAIRE DANS UNE IMAGE DE CIEL
Domaine de l’invention L’invention a trait au domaine du traitement d’image en vue de déterminer la position du soleil dans une image du ciel. L’invention trouve particulièrement application dans la détermination du déplacement de nuages dans un ciel. Le déplacement, caractérisé par la direction et la vitesse, ainsi déterminé peut notamment être utilisée pour l’estimation de l’ensoleillement et de la production d’électricité d’un parc photovoltaïque.
Etat de la technique
La nature irrégulière et intermittente des énergies renouvelables (solaire, hydraulique, etc...) implique des contraintes de production importantes. Ceci explique qu’aujourd’hui encore la production classique d’électricité au moyen de ressources non renouvelables (charbon, nucléaire, etc...), mais totalement contrôlables, est utilisé pour palier la chute de production par énergie renouvelable afin de garantir aux consommateurs un accès garanti à de l’électricité. En outre, ces contraintes et leurs coûts associés sont d’autant plus importants que la chute de production n’est pas anticipée. C’est pourquoi dans le domaine de la production photovoltaïque, il a toujours été porté une attention particulière à la prévision d’ensoleillement, et donc de production photovoltaïque qui en découle, et ce à toute échelle de temps (variabilité annuelle, saisonnière, journalière et infra-journalière).
Par exemple, l’utilisation conjointe de modèles météorologiques et d’images satellite permet de caractériser et prévoir l’évolution de l’ensoleillement aux horizons temporels allant de +lh à +48h. L’advection d’images satellite permet notamment de prévoir la nébulosité aux échéances les plus courtes (jusqu’à +5h), tandis que l’adaptation statistique des modèles météorologiques permet de calculer des prévisions aux échéances plus lointaines.
Cependant, à l’échelle infra-journalière, la production d’énergie sur les parcs photovoltaïques peut chuter brutalement en quelques dizaines de secondes lors du passage de nuages au-dessus du site d’implantation. La reprise est tout aussi soudaine. Cette intermittence est problématique quand l’énergie du parc est utilisée en mode stockage/déstockage en association avec des batteries ou d’autres moyens de production d’énergie renouvelable (centrales hydroélectriques). Les mécanismes d’asservissement sur ces temps de réactions relativement courts peuvent cependant être optimisés si la production d’énergie par les panneaux photovoltaïques est anticipée de quelques minutes à quelques heures. En outre, le couvert nuageux a non seulement une influence sur l’amplitude mais également sur la composition du rayonnement global. La présence de nuages modifie sensiblement l’importance relative des rayonnements diffus et direct. Alors que par ciel clair le rayonnement diffus représente aux mieux 10 à 15% du rayonnement global, les nuages constituent un facteur de diffusion supplémentaire. Par ciel couvert de nuages élevés, l’essentiel du rayonnement global est attribuable au rayonnement diffus.
Or, les techniques à base de modèle météorologique et d’images satellite ne permettant de traiter ni les phénomènes locaux, ni les phénomènes à courte échelle de temps, ni de tenir compte du couvert nuageux. C’est pourquoi l’approche le plus souvent adoptée pour appréhender et anticiper ces phénomènes au voisinage immédiat d’un parc photovoltaïque consiste à acquérir des images du ciel au-dessus du parc depuis le sol, à les traiter pour déterminer la course des nuages dans le ciel et tenter de prévoir sur cette base l’ensoleillement au-dessus du parc photovoltaïque.
Toutefois, cette approche pose des problèmes tant au niveau du matériel mis en œuvre que du traitement d’images réalisé sur le flux d’images acquises. Par exemple acquérir l’image du ciel en une seule prise nécessite l’emploi d’une optique ayant un très grand angle de champ, par exemple de type « fish-eye » qui permet d’obtenir des images hémisphériques englobant la totalité du ciel visible depuis la position de prise de vue. Cependant, la capacité à prendre des photos de très grand angle s’accompagne de nombreux défauts optiques : des distorsions géométriques importantes, des aberrations chromatiques ainsi que d’autres défauts optiques qui mettent ultérieurement en difficulté le traitement d’image et/ou faussent son résultat. Outre la distorsion plus difficile à corriger que pour des optiques classiques, l’un des inconvénients majeurs de ce type d’optique est que la résolution de l’image est très mauvaise lorsque que l’on se rapproche de l’horizon. C’est pourquoi des optiques sophistiquées ont été conçue (e.g. des systèmes catadioptriques panoramique à base de miroir hémisphérique, des miroirs convexes à bande tournante, etc.). Outre leur complexité, ce type d’optique impose l’emploi de pièces mécaniques mobiles (e.g. d’un traqueur pour suivre la course du soleil ou d’un moteur pour entraîner en rotation des éléments optiques), ce qui renchérit le coût global du dispositif, et pose de plus des problèmes de fiabilité et de robustesse.
Concernant le traitement d’image, les techniques les plus répandues pour suivre (« tracking ») des objets en mouvement dans une scène sont de type « Particle Image Velocity » (« PIV »). Cette technique utilise une recherche de ressemblance entre deux images espacées temporellement d’un court laps de temps à partir desquelles il est extrait un champ de vecteurs locaux de déplacement instantané (i.e : direction de mouvement et amplitude de la vitesse). Si cette technique présente des avantages, notamment la capacité à distinguer les zones uniformes des zones fortement variables, elle présente cependant des désavantages nombreux. Notamment, une partie significative de l’image est très souvent non renseignée sur les vecteurs de déplacements et/ou est faussement identifiée comme exempte d’objet en déplacement. Ceci nécessite une calibration très fine du traitement et/ou une post-analyse complexe et coûteuse en paramètres seuils, et donc particulièrement complexe à calibrer. De plus une identification correcte des déplacements par cette technique suppose des déformations faibles des objets (nuages) d’une image à l’autre, ce qui impose des prises d’images à intervalle de temps rapproché. Or, les nuages sont par nature de forme très changeante, notamment parce qu’ils se déforment, parce que de nouveaux nuages apparaissent ou disparaissent, etc... De plus quand les nuages approchent du disque solaire (cas particulièrement intéressant puisqu’en interceptant le disque solaire ces nuages affectent la radiation solaire perçue au sol), leurs propriétés colorimétriques sont souvent très fortement modifiées. Ils deviennent notamment plus clairs dans l’image et mettent donc en difficulté la recherche de ressemblance. Au final, ce type de technique est donc peu robuste pour la détermination de la course des nuages et nécessite d’importantes ressources informatiques pour la mise en œuvre du traitement.
En particulier, le disque solaire dans une image est un élément perturbateur pouvant mettre en difficulté le traitement ultérieur visant à déterminer le déplacement des nuages. En outre, la connaissance de la position du disque solaire dans un flux d’images est un élément important pour prévoir l’ensoleillement (chute ou accroissement). La connaissance précise de cette position nécessite cependant une connaissance précise de la position du dispositif d’acquisition d’image. Cette position est toutefois renseignée manuellement par un opérateur et correspond à la position souhaitée pour ce dernier, et non pas la position réelle du dispositif.
Expose de l’invention
Le but de la présente invention est de proposer un procédé d’estimation de la position du disque solaire dans une image numérique qui permette une calibration automatique et précise de la position de l’appareil d’acquisition d’image. A cet effet, l’invention a pour objet un procédé d’estimation de la position du disque solaire dans une image numérique de ciel acquise au moyen d’un dispositif d’acquisition d’images comprenant une optique apte à focaliser une scène observée sur un plan focal et un capteur numérique disposé dans le plan focal de l’optique, l’estimation de la position du disque solaire étant réalisée en fonction d’un modèle astronomique de l’azimut et de l’élévation du soleil dans le ciel ayant pour paramètre une position prédéterminée (θ[,02,Azr) du dispositif d’imagerie, définie par la position du capteur numérique par rapport à un plan horizontal et par la position d’un axe du capteur par rapport à un azimut (Azr).
Selon l’invention, le procédé comprend une étape de détermination de la position prédéterminée (θ[,θ%,Azr) du dispositif d’imagerie, comportant : - l’acquisition d’un flux d’images numériques de calibration horodatées {IN£} comprenant chacune le disque solaire ; - pour chaque image ΙΝξ dudit flux {IN£} , le calcul d’une position P£ du centre du disque solaire dans l’image IN£ en analysant l’image ΙΝξ ; - le calcul de la position du dispositif d’imagerie (0[, >Azr) comme solution d’un problème d’optimisation selon la relation :
où Ω est un ensemble de recherche, et || || est une fonction de coût quantifiant un écart entre deux vecteurs de R3.
Selon un mode de réalisation, les images sont en couleur, et en ce que le calcul d’une position P£ du centre du disque solaire dans l’image IN£ comprend : - le calcul d’une matrice de luminosité ILcn d’une matrice IN^'R~B = |R(INfi) — B(IN£)\ égale à la valeur absolue de la différence entre des canaux bleu et rouge de l’image IN£ ; et - le calcul de la position P£ dans la matrice ILcn.
En particulier, le calcul d’une position P£ du centre du disque solaire dans l’image IN£ comprend : - le calcul de zones de pixels de la matrice IN^'R~Bcandidats au disque solaire, un pixel (a, b) étant jugé candidat si sa valeur Z£(a,è) dans la matrice ILcn est supérieure à un premier seuil prédéterminé et si sa valeur in^R~B(a,b) dans la matrice INn' est inférieure à un second seuil prédéterminé ; et - la sélection d’une zone candidate comme appartenant au disque solaire en fonction de sa proximité à une position prédéterminée.
Selon un mode de réalisation, le procédé comprend l’estimation de la direction moyenne de déplacements d’objets en mouvement dans un flux d’images numériques {INn} d’une scène comprenant du ciel, les images numériques dudit flux étant munies d’un axe de référence Are^ l’estimation comprenant : - l’acquisition du flux d’images numériques {INn} ; - pour chaque image INn du flux d’images numériques {INn} o le calcul d’une image segmentée IN% identifiant les objets en mouvement dans l’image INn, o le calcul d’un ensemble de zones d’image {Zfcn} respectivement associées à des directions d’analyse {Dk} formant des angles {ak} différents prédéterminés avec l’axe de référence Aref : o pour chaque zone Zkn dudit ensemble de zones {Zkn} , le calcul d’une distribution DDkn des valeurs des pixels de ladite image segmentée IN£ compris dans ladite zone Zkn, la distribution DDkn étant calculée le long d’une direction normale à la direction d’analyse Dk de ladite zone Zkn; - pour chaque direction d’analyse Dk, le calcul d’un critère de variabilité CVkn de distributions {DDkn_p,DDkn_p+1,, DD^n^, DDk>n} calculées pour des zones Zkn associées à ladite direction Dk; et la détermination de la direction moyenne Dkoptn de déplacement des objets dans le flux d’images numériques {INn} comme étant égale à la direction d’analyse pour laquelle le critère de variabilité CVk n est optimum.
En d’autres termes, les inventeurs ont découvert que la zone orientée qui présente la variation temporelle minimale de sa distribution est celle dont la direction correspond à la direction moyenne de déplacement des objets dans une scène. Ce procédé est ainsi bien adapté à la course des nuages dans le ciel. 11 est par ailleurs robuste car il filtre naturellement les phénomènes de haute dynamique/courte durée qui perturbent habituellement les traitements de l’état de la technique. En outre, hormis éventuellement la segmentation des images, les calculs mis en œuvre sont basiques (calcul de distributions et calculs des variations de distribution), et ne nécessitent donc pas de grandes ressources informatiques. De plus, si la segmentation peut nécessiter des ressources ou une calibration, le procédé selon l’invention s’avère cependant très robuste à celle-ci de sorte qu’une segmentation peu précise des images est suffisante. Enfin, il a été constaté que le procédé est robuste également au processus optique de formation des images et renvoie des résultats très satisfaisants, même sur des images réalisées à partir d’un fisheye, y compris lorsqu’on leur applique un redressement sans traitement des aberrations liées à ce type d’optique.
Selon un mode de réalisation, l’acquisition du flux d’image numériques {/n} comprend : l’acquisition de la scène par une dispositif de prise de vue ayant un angle de champ supérieur ou égal à 100° selon la diagonale d’un capteur au format 35mm ou en équivalent au format 35mm, en particulier un angle de champ égal à 180° ; - et préalablement au calcul des distributions DDkn, le redressement des images les images numériques {INn} L’invention exploite ainsi avantageusement les optiques à très angle de champ, par exemple les fisheyes qui permettent en une seule prise de vue d’acquérir l’ensemble de la scène visible, et donc la portion visible du ciel.
Selon un mode de réalisation, la scène est acquise à l’aide d’un dispositif d’imagerie comprenant un système optique apte à former une image de la scène sur un plan focal et un capteur numérique plan disposé sur le plan focal du système optique pour produire le flux d’images numérique {INn} , et le procédé comprend le calcul de l’ensemble de zones d’image {Zkn} pour chaque image INn du flux d’images numériques {INn} , ledit calcul comprenant : le calcul d’un ensemble de rectangles {Æfcn} ayant chacun un petit côté centré sur un même point prédéterminé du plan focal, et ayant leurs grands côtés respectivement parallèles aux directions d’analyse {Dk} ; et le calcul des zones {Zfcn} comme étant égales respectivement à l’intersection des rectangles {Æfcn} avec l’image INn.
En particulier, la scène comprend du ciel, et le point prédéterminé du plan focal est la position du centre du disque solaire sur le plan focal du système optique.
En d’autres termes, le procédé permet de calculer directement le déplacement moyen dans un référentiel d’intérêt, par exemple centré sur l’image du soleil. Il est ainsi déterminé directement les nuages qui se déplacent vers le soleil et qui sont donc susceptible de masquer ultérieurement ce dernier.
Avantageusement, le petit côté de chaque rectangle Rkn de l’ensemble de rectangles {Æ/c,n} est au moins égal au diamètre du halo du disque solaire. Les images sont ainsi découpées en zones qui tiennent compte directement d’un résultat recherché ultérieurement, notamment l’occultation du disque solaire et de son halo par des nuages et donc la variation de l’ensoleillement.
Selon une variante avantageuse, la position du centre du disque solaire est calculée en fonction d’un modèle astronomique de l’azimut et de l’élévation du soleil dans le ciel définissant la position dudit centre sur le plan focal. En d’autres termes, il est possible de déterminer le déplacement des nuages vers le soleil y compris lorsque ce dernier est absent de l’image, par exemple lors qu’il est caché par un nuage, par un arbre, voire même complètement en dehors de l’image.
Selon un mode de réalisation particulier, le critère de de variabilité CVkn pour la direction Dk est calculé en fonction des distributions {DDkn_p,DDkn_p+1, DDkn] associées aux p > 2 dernières images numériques acquises.
Selon un mode de réalisation, le critère de variabilité CVkn pour la direction Dk est calculé selon la relation :
expression dans laquelle j référence la position d’une image numérique dans le flux d’images numériques {INn} , n et n-p référence respectivement la nieme image et la (n-p)ieme image dans le flux d’images numériques {INn} , avec p strictement inférieur àn,k référençant la kieme direction Dk, et ||L|| est une norme du vecteur V de l’espace M1 , / étant la dimension des distributions DDkn.
Les inventeurs ont en effet noté que le critère CVkn présente un minimum global en fonction de l’indice k et que ce minimum global correspond à la direction moyenne de déplacement. Le critère peut être plus au moins raffiné pour mettre en valeur ce minimum global mais même la simple somme de somme des valeurs de différences deux à deux de distribution permet de détecter ce minimum avec précision.
Selon un mode de réalisation, le procédé comporte le calcul d’un masque d’exclusion des éléments immobiles de la scène, et préalablement au calcul des distributions DDkn, ledit masque est appliqué aux images du flux d’images numériques {INn} et/ou aux images segmentées {/N%} . En d’autres termes, en excluant les objets de la scène qui ne comportent aucune information signifiante (arbres, reliefs, bâtiments, etc. . .), le traitement se concentre uniquement sur le ciel et voit donc sa précision renforcée.
Selon un mode de réalisation, le procédé comprend en outre : le calcul d’une vitesse moyenne Vn des objets en mouvement dans le flux d’images numériques {INn} selon ladite direction moyenne de déplacement calculée ; et la prédiction de la position d’un objet en mouvement en fonction de ladite vitesse moyenne Vn et de ladite direction moyenne Dkoptn.
Notamment, les objets en mouvement sont des nuages, et le procédé comprend la prévision d’ensoleillement en fonction de nuages dans le flux d’images numériques {INn} ayant ladite vitesse et ladite direction moyenne. En particulier, le procédé comprend la prévision de l’ensoleillement au-dessus d’un parc photovoltaïque, en particulier la chute ou l’accroissement d’ensoleillement.
En d’autres termes, le procédé exploite avantageusement le fait que la direction moyenne de déplacement est déjà connue pour calculer la vitesse moyenne selon cette direction. Disposant d’une information importante, le calcul de la vitesse est donc facilité, tout type de calcul peut être utilisé à cette fin. Connaissant la direction et la vitesse moyenne, par exemple des nuages dans le ciel, il devient alors possible de faire des prévisions infra-journalières de l’ensoleillement, et donc de la production d’électricité par un parc photovoltaïque.
Selon un mode de réalisation, le calcul de ladite vitesse moyenne Vn comprend : pour chaque image segmentée IN%, le calcul d’une distribution DVk n des valeurs des pixels de l’image segmentée IN£ appartenant à la zone Zk n associée à la direction moyenne déterminée DkoptTl, la distribution DVk n étant calculée le long de la direction moyenne Dkopt Tl déterminée ; le calcul de l’inclinaison d’au moins une structure affine comprise dans une matrice DVn = [DVkoptn_qDVkoptn_q+1... DVk()ptn_iDVkopt1^ issue de la concaténation des distributions DVk n des images entre la (n-q)ieme image du flux jusqu’à la nieme du flux d’images numériques {INn} , où q est un entier supérieur ou égal à 1; et le calcul de la vitesse moyenne Vn en fonction de l’inclinaison calculée.
En d’autres termes, les inventeurs ont constaté qu’il apparaît naturellement des structures affines dans la matrice DVn = [DVkoptin_qDVkoptin-q+1 -DVkoptin-if>Vfcopt,n]· Ces structures affines et parallèles correspondent au déplacement moyen des objets dans le flux d’image et leur pente correspond à la vitesse moyenne des objets selon la direction moyenne de déplacement. Le calcul de la vitesse moyenne est ainsi robuste et se limite à la détermination d’une pente.
Selon un mode de réalisation, le calcul de l’inclinaison de la au moins une structure affine comprise dans une matrice DVn comprend : - l’application de différentes rotations d’angles βη prédéterminés à la matrice DVn ; le calcul d’un histogramme Histnu pour chaque matrice DVnu obtenue par rotation et le calcul du maximum Mn u de cet histogramme ; la détermination d’inclinaison comme étant égal à l’angle βηορ(,η correspondant au maximum des maxima Mn u calculés.
Ce calcul permet de déterminer de manière précise l’inclinaison à l’aide de calcul algébriques simples, et donc peu consommateurs de ressources informatiques. L’invention a également pour objet un système d’estimation de la position du disque solaire dans une image numérique : - une caméra pour l’acquisition de l’image numérique ; - une unité de traitement d’information configurée pour estimer la position du disque solaire dans l’image en fonction d’un modèle astronomique de l’azimut et de l’élévation du soleil dans le ciel ayant pour paramètre une position prédéterminée (0[,0J>Azr) du dispositif d’imagerie, définie par la position (0[, 0£) du capteur numérique par rapport à un plan horizontal et par la position d’un axe du capteur par rapport à un azimut (Azr) ; et ladite unité étant en outre configurée pour déterminer la position prédéterminée (0[, 6r2,Azr) du dispositif d’imagerie : pour chaque image ΙΝξ d’un flux d’images numériques de calibration horodatées {IN£} comprenant chacune le disque solaire {IN£} , en calcul d’une position P£ du centre du disque solaire dans l’image IN£ en analysant l’image IN£ ; en calculant la position du dispositif d’imagerie (0j\0j, Azr) comme solution d’un problème d’optimisation selon la relation :
où Ω est un ensemble de recherche, et || || est une fonction de coût quantifiant un écart entre deux vecteurs de M3.
En particulier, le système est configuré pour mettre en œuvre un procédé du type précité. L’invention a également pour objet un produit programme d’ordinateur enregistré sur support informatique utilisable par un ordinateur comprenant des instructions pour l’exécution d’un procédé d’estimation de la position du disque solaire dans une image numérique de ciel acquise au moyen d’un dispositif d’acquisition d’images comprenant une optique apte à focaliser une scène observée sur un plan focal et un capteur numérique disposé dans le plan focal de l’optique, l’estimation de la position du disque solaire étant réalisée en fonction d’un modèle astronomique de l’azimut et de l’élévation du soleil dans le ciel ayant pour paramètre une position prédéterminée (θ[, θζ,Αζ1") du dispositif d’imagerie, définie par la position (θ[,θ£) du capteur numérique par rapport à un plan horizontal et par la position d’un axe du capteur par rapport à un azimut (Azr), le procédé comprenant une étape de détermination de la position prédéterminée (β^,θ^,Azr) du dispositif d’imagerie, comportant : l’acquisition d’un flux d’images numériques de calibration horodatées {ΙΝξ} comprenant chacune le disque solaire ; pour chaque image IN£ dudit flux {IN£} , le calcul d’une position P£ du centre du disque solaire dans l’image IN£ en analysant l’image ΙΝξ ; - le calcul de la position du dispositif d’imagerie (0[, >Azr) comme solution d’un problème d’optimisation selon la relation :
où Ω est un ensemble de recherche, et || || est une fonction de coût quantifiant un écart entre deux vecteurs de M3.
En particulier, le support comprend des instructions pour la mise en œuvre d’un procédé du type précité.
BREVE DESCRIPTION DES FIGURES L’invention sera mieux comprise à la lecture de la description qui va suivre, donnée uniquement à titre d’exemple, et faite en relation avec les dessins annexés dans lesquels : - la figure 1 est une vue schématique d’un système de prédiction d’ensoleillement selon l’invention ; les figures 2A-2D sont des organigrammes d’un procédé de prévision d’ensoleillement prédiction de la production photovoltaïque mis en œuvre, par exemple mis en œuvre par le système de la figure 1 ; la figure 3 est une image, transformée en noir et blanc, d’une image acquise par le système de la figure 1 ; - la figure 4 est une image issue d’une segmentation, e.g. de type Otsu, de la figure 3 pour la détermination d’objets non météorologiques ; la figure 5 est un profil illustrant la localisation d’une portion de ciel dans l’image de la figure 3, issu de la segmentation de la figure 4 ; la figure 6 est une image numérique d’un masque d’exclusion obtenu à partir du profil de la figure 5 ; la figure 7 est une image numérique d’un masque d’exclusion obtenu par application d’un masque de dilation circulaire à l’image de la figure 6 ; - la figure 8 est une image numérique illustrant l’application du masque de la figure 7 à l’image numérique de la figure 3 ; - la figure 9 est une image numérique illustrant l’application d’un masque comprenant en outre le disque solaire à l’image numérique de la figure 3 ; - la figure 10 est une image numérique illustrant une segmentation pour la détection des nuages dans l’image de la figure 9 ; - la figure 11 est une image illustrant un rectangle de détection selon l’invention ; - les figures 12A-C sont des figures illustrant différent type de zones de détection selon l’invention ; - les figures 13A-E sont des images illustrant des flux de distributions de déplacement pour différents angles de détection ; - la figure 14 est une vue schématique illustrant graphiquement le calcul d’un flux de distribution ; - les figures 15 à 17 sont des tracées illustrant différents critères de variabilité de flux en fonction de l’angle de détection ; la figure 18 illustre une comparaison de la détection automatique de la direction moyenne de déplacement selon l’invention et une détection manuelle de cette direction ; la figure 19 est une vue schématique illustrant graphiquement le calcul d’une distribution de vitesse selon l’invention ; la figure 20 est une image illustrant graphiquement des structures affines d’un flux de distributions de vitesse ; la figure 21 est une image illustrant une rotation des structures des figures de la figure 20 ; la figure 22 est un tracé d’un critère d’une distribution de vitesse en fonction d’un angle de rotation de cette distribution; la figure 23 illustre la comparaison entre la détermination automatique de l’angle de rotation correspondant à la vitesse de déplacement des nuages et la détermination manuelle de cet angle ; la figure 24 est une image illustrant une partie du flux de distribution de vitesse ; les figures 25A-G sont des tracés illustrant des échéances temporelles d’arrivée de nuage sur le disque solaire pour différentes distances en pixel du soleil, les figures 26 et 27 sont des tracées superposant les instants estimés d’arrivée des nuages sur le soleil, la figure 27 illustrant en outre des mesures d’éclairement par un pyranomètre ; et les figures 28 et 29 sont des images illustrant le procédé de calibration de la position du dispositif de prise de vue selon l’invention.
Description detaillee de l’invention A. Notations
Dans la suite du document, les conventions mathématiques classiques de l’algèbre linéaire sont appliquées pour des matrices réelles de M/1xS, où A et B sont des entiers positifs, avec notamment pour des matrices M et N de RAxB : M - (ma,ù)ae[1A]jùe[1B] \M\ - (\ma,b\) ae[lA],be[l,B]
La notation Μ ® N se rapporte à l’application d’un masque d’exclusion Ma la matrice N. Plus particulièrement, l’application d’un traitement sur l’objet M ® N signifie que seuls les pixels de la matrice N correspondant à des pixels du masque M ayant une première valeur sont pris en compte pour le traitement. Par exemple ma b = 1 si le pixel est pris en compte ou 0 s’il est exclu du traitement. De même, lorsqu’il est fait référence à l’objet Μ ® N seul (la « matrice M ® JV »), il est fait référence à l’ensemble de pixels de N qui ne sont pas exclus par le masque M La manière de prendre ou non en compte un pixel de la matrice N en fonction du masque M peut être réalisée de nombreuses manières connues en soit, la notion de masque étant classique en traitement d’images. Par convention, dans la demande la « matrice » M ® N est représentée graphiquement par l’image de N modifiée en réglant ses pixels en noir lorsqu’ils correspondent aux pixels du masque M de valeur 1. B. MATERIEL·
En se référant à la figure 1, un dispositif 10 d’imagerie et de traitement d’image selon l’invention, directement implanté sur un parc photovoltaïque, comporte un appareil photographique ou une caméra 12, par exemple de type reflex, muni d’un capteur numérique CMOS ou CCD 14 pour l’acquisition d’image dans le rayonnement visible, par exemple au format 35mm, et d’un objectif 16 à très grand angle de champ supérieur à 100° sur la diagonale du capteur 14, avantageusement de 180° (e.g. fisheye). L’objectif produit ainsi sur le capteur 14 une image optique circulaire de l’espace devant l’objectif 16 et donc une image comprenant l’ensemble du ciel visible depuis la position du dispositif 10. Le capteur 14 est agencé sur le plan focal 18 de l’objectif 16 et centré sur l’axe optique 20 de ce dernier. L’appareil 12 comprend par ailleurs les circuits 22 nécessaires au pilotage de l’acquisition des images et le prétraitement des données issues du capteur (e.g. dématriçage) pour produire une image à partir des signaux générés depuis les photosites sensibles du capteur 14, usuellement un processeur « DSP » (pour « digital signal processor», ainsi qu’une horloge permettant d’horodater les images acquises. L’acquisition des images est réalisée avec une période d’acquisition Te, l’appareil comportant à cet effet un intervallomètre interne ou le dispositif comportant un intervallomètre externe d’une manière connu en soi. Par exemple, l’appareil 12 est un appareil Canon de référence 6D équipé d’un objectif Canon EF 8-15 mm f/4 L réglé sur la focale 8mm, piloté par un intervallomètre externe de référence DigiSnap 2700 de la société Harbonics.
Le dispositif 10 comporte par ailleurs une unité de traitement d’informations 24 connectée à l’appareil 12, par un câble ou par une communication sans fil, pour recevoir les images acquises par l’appareil 12, par exemple au format brut (après le dématriçage mais avant tout traitement informatique, correspondant au format « RAW » chez Canon) ou dans un format de lecture (e.g. jpeg), et pour leur appliquer un traitement mis en œuvre par ordinateur selon l’invention. L’unité 24 est par exemple un ordinateur personnel, une tablette, un smartphone, ou plus généralement tout système à base de microprocesseur s), notamment de type DSP (« digital signal processor »), à base de circuits de type FPGA, à base de circuits mixant ces types de technologie, etc., apte à mettre en œuvre l’invention. L’unité 24 est notamment pourvue de l’ensemble des mémoires (RAM, ROM, cache, mémoire de masse,...) pour la mémorisation des instructions mettant en œuvre le procédé selon l’invention, des paramètres et réglages utiles à cette mise en œuvre et pour la mémorisation des résultats des calculs intermédiaires et finaux.
Optionnellement, l’appareil photographique 12 est avantageusement logé dans un boîtier de protection 26 muni d’une fenêtre transparente 28 en face de l’objectif 16 pour le protéger des intempéries et ainsi permettre de laisser le dispositif à demeure sur le site choisi. Dans une variante, l’unité 24 est également logée dans le boîtier de protection 26 et communique, via par exemple un réseau sans fil (téléphonique, satellitaire, etc...), le résultat du procédé qu’elle met en œuvre à un site distant où la production du parc photovoltaïque est gérée.
En fonctionnement, un flux d’images numériques horodatées {INn} est donc produit par l’appareil 12, une image numérique INn étant acquise par exemple toutes les Te secondes, l’entier n représentant le nlème instant d’acquisition nxTe depuis un instant initial n = 1, et donc la //'eme image dans le flux d’images {INn} . En variante, l’acquisition est réalisée de manière non régulière avec un pas d’acquisition non constant dans le temps.
Chaque image INn comprend A lignes de pixels par B colonnes de pixels et comprend trois canaux d’un espace colorimétrique, par exemple les canaux rouge (r), vert (v) et bleu (b) codés respectivement par trois matrices R, F, B de RAxB. Les valeurs numériques composant les matrices sont codées sur un nombre bits prédéterminé, par exemple 8 bits, correspondant à 256 niveaux de valeurs (de 0 à 255). Par la suite, les canaux r,v,b d’une image en couleur quelconque « IN » sont respectivement notés R(IN), V(IN) et B(IN).
Dans ce qui suit l’espace colorimétrique considérée sera l’espace RVB, étant entendu que l’invention s’applique à d’autres espaces colorimétriques, par exemple l’espace HSV (teinte, saturation, intensité lumineuse). De même, il est décrit un nombre de bit de codage égal à 8, correspondant à des valeurs comprises entre 0 et 255. Bien entendu, l’invention s’applique à un nombre de bits quelconque. C. Traitement informatique
Il va à présent être décrit, en relation avec les organigrammes des figures 2A-2E, un procédé selon un mode de réalisation de prévision de la production d’électricité d’un parc photovoltaïque, et plus particulièrement la chute d’ensoleillement sur le parc photovoltaïque entraînant une chute de la production, mis en œuvre par le dispositif de la figure 1, ce procédé se fondant sur l’estimation selon l’invention d’une direction moyenne de déplacement et d’une vitesse moyenne des nuages dans le ciel. C.l. Calibration
Le procédé débute par une phase 40 de positionnement et de calibration du dispositif d’imagerie 10.
Une première étape 400 consiste à charger dans la mémoire de l’unité 24 l’ensemble des instructions et paramètres utilisés par le procédé, notamment : - un modèle astronomique de la position du soleil dans le ciel, par exemple déterminant l’azimut et l’élévation du soleil dans le ciel en fonction de la longitude, la latitude, de l’heure et du jour de l’année. Ce modèle est par exemple le « solar position calculator » de la NOAA (« National Oceanic and Atmospheric Administration ») ; - un modèle de la formation d’image sur le capteur 14 par l’objectif 16. Pour le fisheye, il s’agit par exemple du modèle décrit dans le document « Géométrie modelling and calibration of fisheye lens caméra System », Schwalbe E., In Proc. 2nd Panaramic Photogrammetry Worshop, Int. Archives of Photogrammetry and Remote Sensing (vol. 36, No. Part5, p. W8), février 2005. Ce modèle permet notamment de modéliser l’image d’un objet photographié par le dispositif 12 sur le plan focal de l’objectif 16, et donc sur le capteur 14.
En connaissant la position θ2,Αζ) du capteur 14 dans un repère terrestre, et donc la position du capteur par rapport au soleil, la date et l’heure, la combinaison du modèle astronomique et du modèle de la formation d’image permet de calculer l’image du centre du soleil sur le plan focal de l’objectif, et donc le centre du disque solaire dans les images produites par l’appareil 12. Cette combinaison est ci-après désignée par les termes « modèle astronomique optique ».
En 402, un opérateur positionne l’appareil 12 sur un site choisi avec l’objectif visant le ciel. La position de l’appareil est ainsi défini par trois paramètres (01# θ2,Αζ), le couple d’angles (θ±, θ2) définissant la position du capteur 14 par rapport à un plan horizontal, et le paramètre Az représentant l’angle que forme l’axe passant par les centres des bords haut et bas du capteur avec un azimut, par exemple l’azimut nord. L’opérateur cherche notamment à positionner le capteur 14 sur un plan horizontal avec l’azimut nord correspondant au milieu du bord supérieur des images produites par l’appareil 12. La position souhaitée, i.e. idéale, de l’appareil est alors égale à (θνθ2,Αζ) = (0,0,0). Les coordonnées géographiques, e.g. la longitude et la latitude, du dispositif sont par ailleurs mémorisées dans l’unité 24.
Evidemment, il est recherché à positionner le dispositif dans une zone dégagée de manière à maximiser la portion de ciel dans les images. Toutefois, le dispositif peut être implanté dans une zone très contrainte en termes d’éléments non pertinents présents dans le champ de l’appareil (ou « objet non météorologiques », ou « onm », par exemple des arbres, des constructions, des pilonnes, etc...), l’invention permettant en effet une prévision robuste y compris sur des portions minimes de ciel et y compris avec le soleil caché.
Une fois le dispositif 10 positionné, une calibration 404 est mise en œuvre. Notamment, le positionnement parfait de l’appareil est difficile à obtenir, de sorte qu’il présente une position réelle (θ[, 02,Azr) différente de la position théorique (θ1,θ2,Αζ). La position du dispositif étant utilisée lors d’étapes de traitement selon l’invention, une correction de l’écart entre les positions réelle et théorique est mise en œuvre. L’étape de calibration 402 a pour objectif de déterminer les paramètres de cette correction et est décrite plus en détail par la suite. La position corrigée (θ[, 02,Azr) est mémorisée par l’unité 24 et est celle utilisée par le modèle astronomique optique pour calculer le centre du disque solaire dans les images produites par l’appareil 12.
La phase initiale 40 se poursuit par le calcul, en 406, d’un masque d’exclusion des objets non météorologiques de la scène présents dans le flux d’images acquis. Ce masque permet de définir les pixels des images sur lesquels portera le traitement pour la prévision de l’ensoleillement. Etant dans une phase d’initialisation et de calibration, il est possible de mettre en œuvre toute méthode connue permettant de déterminer ce type de masque, y compris les méthodes les plus sophistiquées et consommatrices de temps et de ressources matérielles, voire même avec l’aide de l’opérateur qui désigne à la main tout ou partie de ces objets.
De manière avantageuse, la détermination du masque d’exclusion se fonde sur les spécificités des images acquises à savoir celle du ciel, ce qui permet une détermination automatique et rapide du masque. Notamment, en journée, le ciel « clair » (portion de ciel sans nuage et sans le disque solaire, i.e. les portions de ciel bleu), les nuages, le disque solaire et les objets non météorologiques ont des luminosités différentes, comme illustré par exemple à la figure 3.
La détermination du masque d’exclusion comporte ainsi : a) l’acquisition par l’appareil 12 d’une image numérique INm en journée, par exemple par temps ensoleillé ou légèrement nuageux avec le soleil haut dans le ciel, e.g. à son zénith ; b) le calcul d’une matrice de luminosité ILm de l’image INm, d’une manière connue en soi, par exemple selon la formule ILm = - X (R + V + B), obtenant ainsi une image en niveau de gris telle qu’illustrée à la figure 3 ; c) la détermination d’un premier masque Mx G RAxB correspondant à l’image optique de la scène sur le capteur, à savoir m1(a, b) = 1 si le pixel (a, b) de INm appartient à l’image de la scène sur le capteur 14 (i.e. au disque fisheye sur l’image) et 0 sinon ; d) l’application d’un algorithme de segmentation d’histogramme de luminosité, tel que par exemple décrit dans le document « A threshold sélection method form gray-level histograms», de Otsu N, Automatica, 11(285-296), 23-27, 1975. Ce type d’algorithme consiste à diviser l’histogramme de la matrice Mt 0 ILm en une pluralité de parties homogènes du point de vue de leurs valeurs, par exemple quatre parties en calculant trois seuils SI, S2, S3 tels que: (51,52,53) = argmin (ff(0,sl) + a(sl + 1, s2) + a(s2 + l,s3) + a(s3 + 1,255)) (sl,s2,S3)e[0,255]3 où σ(η1,η2) est l’écart type des valeurs de l’histogramme d’abscisses comprises entre ni et n2. En adoptant 5 niveaux de gris NG1=0 <NG2 < NG3< NG4 NG5= 255, l’image M1 0 résultante de cette segmentation (il^tu{a,b) = N G1 si ilm(a,b) < 51, il^tu(a,b) = NG2 si 51 < ilm < 52, etc.) est illustré à la figure 4. Ce type de problème est simple de résolution à l’aide d’un solveur approprié, par exemple mis en œuvre par l’unité 24 ; e) Pour chacun des angles d’un cercle trigonométrique centré sur le centre de l’image (entre 0° et 360° avec un pas angulaire p prédéterminé, par exemple de 1°), il est déterminé le pixel de la matrice ILm le plus proche du centre du cercle dont la luminosité iZm(a,b) est inférieure au seuil SI qui définit la classe de luminosité la plus sombre issue de la segmentation. Il est ainsi obtenu un profil, par exemple mémorisé dans un tableau de 360lp valeurs d’angles, de la distance au centre de l’image qui est la plus petite distance identifiée d’un pixel sombre, le profil résultant de l’image Mt 0 ILm étant illustré à la figure 5 ; f) la détermination d’un second masque M2 e Ey4xBavec m2(a,b) égal à 0 pour tous les pixels au-delà et sur le profil et égal à 1 sinon, le masque M2 étant illustré à la figure 6.
Le masque M2 est dans une première variante celui utilisé lors du traitement ultérieur de prévision d’ensoleillement. Toutefois à la figure 6, on constate que les limites de la zone masquée (en noir) sont très proches par endroits des cimes des arbres et des contours des branches. On se doute qu’en cas de vent, ce masque ne permettrait pas d’exclure les branches qui se déplaceraient et apparaîtraient alors dans la partie visible (en blanc) du masque. Dans une seconde variante, une correction supplémentaire du masque M2 est appliquée pour tenir compte de ce phénomène. Notamment, un algorithme de dilatation est mis en œuvre sur le masque M2 avec un rayon suffisant pour englober au moins le feuillage, par exemple au moyen d’un masque de dilatation circulaire ayant un rayon compris entre 10 et 50 pixels. Plus particulièrement, pour chacun des pixels préalablement défini comme masqué du masque M2 (m2(a, b) = 0), un masque circulaire est centrée sur ce pixel et tous autres pixels de M2 contenu dans le masque circulaire sont réglés sur 0. Le masque M3 résultant de la dilation du masque M2 à l’aide d’un masque de dilation circulaire de rayon égal par exemple à 10 pixels est illustré à la figure 7. En se référant à la représentation de la matrice M3 0 ILm de la figure 8, on peut ainsi constater que le masque d’exclusion permet d’identifier les parties utiles de l’image (ici le ciel, le nuage et le soleil) sur lesquels porteront les traitements ultérieurs mis en œuvre dans le cadre de la prévision de la chute d’ensoleillement sur le parc photovoltaïque.
Le masque d’exclusion correspondant aux objets non météorologiques Monm = M2 ou Monm = M3 est alors stockée dans une mémoire de l’unité 24.
Une fois la phase initiale de positionnement et calibration terminée, le dispositif est mis en service. L’appareil 12 produit le flux d’images numériques {INn} qui est traité en 50 par l’unité d’information 24 pour l’estimation de l’ensoleillement. Le procédé ayant pour but de déterminer les variations infra-journalières d’ensoleillement, le traitement complet est mis en œuvre avantageusement uniquement lorsqu’il existe du ciel clair. Par exemple, l’unité 24 détermine la luminosité moyenne de chaque image INn et l’unité 24 déclenche le traitement si la luminosité moyenne (e.g. de l’image Monm 0 ILm) dépasse un seuil prédéterminé, signifiant que le ciel n’est pas couvert. D’autres tests, ou aucun test, de déclenchement peuvent bien entendu être mis en œuvre. L’étape de traitement informatique 50 comprend, pour chaque nouvelle image numérique acquise INn, une étape 500 de prétraitement de l’image, suivie d’une étape 502 d’estimation de la direction moyenne des nuages à l’instant d’acquisition n de l’image, suivi d’une étape 504 d’estimation de la vitesse moyenne des nuages audit instant n, suivi par une étape 506 de prévision de la chute de l’ensoleillement au dessus du système 12, et donc du parc où ce dernier est implanté, en fonction de la direction et de la vitesse moyenne calculée à l’instant d’acquisition n. C.2. Prétraitement de l’image lNn
Afin de faciliter la segmentation ultérieure entre ciel clair et nuages, le procédé débute, en 500-2, par l’exclusion supplémentaire du disque solaire si ce dernier est présent au moins partiellement dans l’image INn nouvelle acquise, ce qui ne laisse donc subsister que les portions de l’image acquise INn correspondant au ciel clair et aux nuages. Plus particulièrement, l’unité 24 calcule la position théorique Pn du centre du disque solaire dans le plan focal et donc le cas échéant dans l’image INn en fonction du modèle astronomique optique, puis détermine un disque de rayon rs ayant pour centre la position Pn. Le rayon rs, paramètre enregistré dans l’unité 24, est une valeur déterminée par exemple lors d’une campagne de test antérieure et est supérieure ou égale au plus grand rayon de disque solaire observé au cours de l’année, ou d’une période définie de l’année, sur toute l’année ou sur plusieurs années. D’autres méthodes de calcul dudit rayon sont évidemment possibles. Le rayon rs permet donc à d’exclure de manière certaine l’image du soleil dans l’image numérique INn sans pour autant exclure une zone trop importante. Cette valeur, qui dépend de la résolution de l’image, est par exemple égale à 80 pixels.
Un masque d’exclusion M4n de M/1xB est alors calculé par l’unité 24. Ce masque est égal au masque Monm lorsque le disque solaire n’est pas compris dans l’image. Dans le cas contraire, le nouveau masque d’exclusion M4n correspond à l’exclusion des objets non météorologiques et des pixels du disque solaire compris dans l’image INn, à savoir m4„(a, b) = 1 si m3(a, b) = 1 ou si (a, b) appartient au disque solaire, et m4n(a, b) = 0 sinon. L’image M4n 0 ILn correspond ainsi uniquement aux portions de ciel clair et de nuages de la scène photographiée et est illustrée à la figure 9.
Une segmentation de l’image M4n 0 INn est ensuite mise en œuvre, en 500-4, pour identifier dans celle-ci les portions de ciel clair et les nuages.
Avantageusement, la matrice M4 n 0 INB~B est utilisée pour cette segmentation, où IN*~b = |R(/Nn ) — B(jNn )|. Comme décrit ci-dessous en relation avec la calibration de la position réelle de l’appareil 12, la matrice INB~B, issue de la différence des canaux rouge et bleu, fait apparaître de manière contrastée les portions de ciel clair et les nuages. Un simple seuillage sur la valeur centrale de l’histogramme de cette matrice peut ainsi suffire à segmenter l’image. Un algorithme de type Otsu tel que décrit précédemment, avec un unique seuil, est par exemple alternativement mis en œuvre pour distinguer dans l’histogramme la classe du ciel clair et la classe des nuages. Ceci permet une segmentation plus précise, notamment concernant les portions de nuages. Suite à cette segmentation, une matrice binaire M4 n 0 IN£cn est ainsi calculée (« CCN » pour « ciel clair/nuage ») avec par exemple in^CN (a, b) = 1 si le pixel (a, b) appartient à un nuage et 0 si le pixel (a, b) appartient à une portion de ciel clair. Cette matrice est illustrée à la figure 10.
Le prétraitement 500 se poursuit en 500-6 par la correction de la distorsion géométrique induite par l’objectif 16, notamment la correction du masques M4n et de l’image binaire INfiCL. Ce type de traitement, également connu sous le terme de « redressement », vise à ce qu’une droite dans la scène photographiée soit également une droite dans l’image redressée.
De manière privilégiée, le redressement mis en œuvre se fonde sur le fait que l’invention est peu sensible aux effets de perspectives dans l’image, et utilise donc avantageusement le redressement simplifié de la projection d’une scène 3D sur un plan 2D par un fisheye, par exemple celui décrit dans le document de Schwalbe E. cité précédemment Ce redressement, qui ne corrige que les déformations de l’image liées à la projection au travers de l’objectif, est donc très simple et rapide à calculer.
Dans ce qui suit, les portions de ciel clair n’étant plus utilisées en tant que telles pour la détermination de la direction des nuages, le procédé se poursuit par exemple par la génération d’une image binaire redressée, notée IN£, identifiant les nuages de l’image INn, à savoir in£(a, b) = 1 si le pixel (a, b) appartient à un nuage et 0 sinon. Cette image inclut donc le masque M4 n.
Enfin, l’unité 24 calcule en 500-8 la position Sn = (xn,yn) du centre du disque solaire dans l’image redressée IN£. Le prétraitement mis en œuvre par l’unité 24 calcule ainsi pour chaque image numérique INn acquise par le dispositif 12, une image binaire redressée IN%, ainsi que la position Sn du centre du disque solaire dans l’image redressée.
On notera que les dimensions des images redressées ne sont pas nécessairement égales à celles des images issues du capteur. Pour des raisons de concision, il sera par la suite fait cependant cette hypothèse, les matrices et calculs pouvant bien évidemment être définis pour des dimensions différentes. La même hypothèse est utilisée pour la rotation d’une matrice.
Le procédé continue par l’estimation de la course moyenne des nuages à l’instant n, et plus particulièrement par l’estimation 502 de leur direction moyenne, suivi de l’estimation 504 de leur vitesse moyenne selon la direction moyenne. C.3. Estimation de la direction moyenne des nuages L’estimation de la direction moyenne des nuages consiste dans une première étape à découper la nouvelle image IN£ issue du prétraitement de l’image INn en une pluralité de zones, chacune associée à une direction particulière de recherche Dk prédéterminée et constante dans le flux d’images acquises {INn} . Puis, dans une seconde étape, l’estimation détermine l’évolution des nuages en fonction du temps dans chacune des zones. La direction Dk de la zone de variation nuageuse minimale correspond alors à la direction moyenne des nuages. De manière avantageuse, l’estimation selon l’invention est configurée pour tenir compte de la prévision d’ensoleillement recherchée.
En se référant à la figure 11 illustrant un rectangle sur une image M4n®INn redressée, l’étape 502 d’estimation de la direction moyenne des nuages débute, en 502-2, par la définition d’un ensemble de rectangles de détection {/?fcn} de longueur et de largeur en pixels identiques mais d’orientations différentes. Plus particulièrement, les rectangles de détection Rkn ont : - une largeur en pixels / égale au moins égale au diamètre du halo solaire, noté D%. Le halo solaire, également appelé « cône solaire », correspond à une zone de luminosité très proche de celle du disque solaire et entourant ce dernier, dont l’occultation par des nuages induit donc également une chute de l’ensoleillement. La largeur / est par exemple déterminée lors d’une campagne de test antérieure, mémorisée dans l’unité 24 et est supérieure ou égale au diamètre maximal du halo solaire observé lors de cette campagne. Par exemple, pour une image IN£ de 500 pixels par 500 pixels, la largeur / est égale à 160 pixels ; - un petit côté centré sur le centre Sn du disque solaire. Le rectangle permet ainsi de tenir compte de l’ensemble des nuages susceptibles d’intercepter ultérieurement le halo solaire D%, et ainsi faire varier l’ensoleillement ; - une longueur en pixels L, par exemple égale à la plus grande dimension des images IN£ (e.g. 500 pixels dans les exemples illustrés sur les figures). Comme détaillé ci-après, la longueur des rectangles définit la fenêtre temporelle de prévision des variations de l’ensoleillement. Notamment, plus la longueur L est importante et plus les nuages éloignés du soleil sont pris en compte pour la prévision d’ensoleillement, notamment la chute d’ensoleillement; - une direction prédéfinie Dk, parallèle à la longueur L, repéré par un angle ak par rapport à un axe prédéterminé de l’image, par exemple un axe horizontal des images (e.g. leur bord inférieur). Les directions Dk sont par exemple régulièrement espacées angulairement entre 0° et 360°, avec un pas angulaire Θ qui règle la précision souhaitée sur la direction moyenne recherchée. Par exemple, le pas Θ est choisi égal à 1°, résultant en K = 360 directions Dk fixes étudiées.
Comme illustré aux figures 12A-C, en fonction de la position du disque solaire dans l’image redressée, un rectangle peut être compris entièrement dans l’image (figure 12A), ou partiellement dans l’image si le disque solaire est dans l’image et proche d’un bord de celle-ci (figure 12B) ou si le disque solaire est hors de l’image et proche d’un bord de l’image (figure 12C).
Le procédé consiste, dans une étape suivante 502-4, à déterminer, pour chaque zone Zfc n de l’image binaire redressée IN£, comprise dans un rectangle Rkn, la distribution DDkn des pixels (DD pour « distribution déplacement ») correspondant aux nuages le long de la largeur / du rectangle Rk n, et donc dans une direction normale à la direction Dk associée à ce rectangle.
Par exemple, pour chaque rectangle de détection Rk n, le procédé : produit une image binaire rectangulaire lRkiU correspondant au contenu du rectangle Rkn (i.e. la zone Zfcn, le cas échéant complétée d’une portion de pixels à 0 correspondant à la zone du rectangle Rkn non comprise dans l’image IN£, tels qu’illustré aux figures 12B et 12C) ; applique une rotation d’angle 90° — ak afin d’obtenir une image rectangulaire IRkn de l x L pixels, le petit côté du rectangle étant ainsi horizontal ; calcule la distribution DDkn. Par exemple, la distribution DDkn est un histogramme comptabilisant pour chaque colonne de l’image lRkn le nombre de pixels de la colonne qui sont égaux à 1, i.e. correspondant à des nuages.
Pour chaque image numérique INn acquise par le dispositif 10, l’unité 24 calcule donc un ensemble de distributions {DDkn} respectivement associées à des directions fixes {Dk} . L’unité 24 produit donc, pour le flux d’images {INn} , K flux de distributions, par exemple noté pour le £ieme flux FDDkn = {DDkl, DDk2, —, DDkn_1,DDkn}. Des exemples de flux de distribution FDDkn sont illustrés graphiquement aux figures 13A-13J qui correspondent respectivement aux angles ak égaux à 0°, 10°, 20°, 30°, 40°, 50°, 60°, 70°, 80° et 90°. En se référant à la figure 14, la représentation graphique consiste à transformer l’échelle des valeurs des distributions (de 0 à I concernant des histogrammes) en échelle de gris, ce qui donne une représentation graphique d’une distribution sous la forme d’une étiquette, puis à juxtaposer les étiquettes en fonction des valeurs croissantes de n.
Selon l’invention, l’analyse de la variabilité des flux de distributions FDDkn permet de déterminer la direction moyenne des nuages à l’instant n. A cet effet, dans une étape suivante 502-6 du procédé, la variabilité à l’instant n est calculée en fonction des p dernières distributions calculées {DDkn_p, DDfcn_p+1, ...,DDkn_v DDkn], i.e. sur une fenêtre temporelle [(n — p) xTe;nx Te\, où p est entier prédéterminé. La durée de la fenêtre est par exemple comprise entre 4 secondes et 44 secondes, ce qui correspond à p= 1, 2,.., ou 11 pour un pas d’acquisition Te de 4 secondes. La variabilité est par exemple calculée par tout critère approprié CVkn, notamment un critère de ressemblance ou de dissemblance entre les distributions comprises entre l’instant d’acquisition n et l’instant antérieur d’acquisition n-p.
Avantageusement, le critère CVkn est égal à la somme des ressemblances entre deux distributions espacées de m > 1 pas de temps, à savoir un critère de dissemblance CVkn pour m = 1 selon la relation :
(1)
Les inventeurs ont découvert qu’un tel critère comprend deux optimums (deux maxima si critère de ressemblance ou deux minima si critère de dissemblance) en fonction de l’indice k. Ces angles correspondent respectivement à la direction l’angle des nuages approchant le disque solaire et l’angle des pixels nuages s’éloignant du disque solaire (i.e. correspondant à un angle ak et un angle ak + 180°). Autrement dit si les nuages viennent du nord-ouest, la ressemblance des distributions est maximale pour un angle de détection akn de 135° mais également pour l’angle 135°+180°= 215° (direction sud-est) puisque, selon cette direction, les mêmes nuages qui, une fois passés sous le disque solaire, s’en éloignent avec la même direction.
Le procédé détermine ainsi l’indice k, correspondant à l’angle ak, qui maximise (resp. minimise) le critère CVkn entre 0° et 180°, ainsi que l’indice k\ correspondant à l’angle ak + 180°, qui maximise (resp. minimise) le critère CVkn entre 180° et 360°. Par exemple, pour un pas d’angle Θ de 1°, le procédé détermine l’indice k=kopt compris entre 0 et 180 qui maximise (resp. minimise) le critère CVkn = CVkn + CVk+18Q n. La direction Dkoptn correspondant à l’indice kopt déterminé par l’unité 24 pour l’instant n, est ainsi la direction moyenne de déplacement des nuages à l’instant n.
La figure 15 illustre des critères CVkn se fondant sur la relation (1) pour différentes écartements temporels entre les distributions. Les ronds indiquent la valeur de l’angle (Xkopt Pour laquelle le critère CVkn est minimal, et donc la direction Dkoptn de propagation des nuages à l’instant n finalement identifié par le procédé selon l’invention. On constate que le minimum de ces sommes est faiblement sensible à l’écart temporel entre deux distributions. Notamment, l’angle de propagation des nuages varie entre 30 et 28 degrés quel que soit la durée de la fenêtre. D’autres critères CVkn que celui de la relation (1) peuvent bien entendu être utilisés. Notamment d’autres types de norme que la norme L1 d’un vecteur d’un vecteur peuvent être utilisées. De même, en variante un critère R2 ou de Nash est utilisé. Les résultats de ces critères de ressemblances sont illustrés respectivement aux figures 16 et 17. Tous deux sont d’autant meilleurs qu’ils approchent la valeur de 1, le procédé recherchant cette fois un maximum et non un minimum comme précédemment. On constate que ces deux critères ont un maximum peu marqué, ce qui peut rendre sa détermination plus difficile. Le critère selon la relation (1), qui présente un minimum très marqué, est donc privilégié.
La figure 18 illustre une direction de propagation déterminée manuellement par un opérateur. Sur deux images espacées de 7 pas de temps (soit 28 secondes) l’opérateur a dessiné à la main trois vecteurs de déplacements de nuages. L’opérateur a identifié dans les deux images deux zones correspondant aux mêmes nuages. L’origine d’un vecteur correspond ainsi à un pixel de nuage dans la première image et l’extrémité du vecteur correspond au pixel correspondant dans la seconde image. Ces trois mesures de direction de déplacement des nuages sont illustrées par les trois rapporteurs et donnent des valeurs sensiblement identiques et cohérentes avec la direction déterminée par le procédé selon l’invention (30°), étant entendu que le positionnement des vecteurs par l’opérateur est incertain à quelques pixels près en raison de la déformation des nuages.
La direction moyenne de propagation des nuages Dkoptn étant ainsi connue à chaque instant n, il est possible de mettre en œuvre une détermination de leur vitesse moyenne, selon tout traitement connu de l’état de la technique. Un calcul selon l’invention de cette vitesse est à présent décrit. C.4. Estimation de la vitesse moyenne de déplacement des nuages
Une fois déterminée la direction moyenne de déplacement Dkoptn des nuages à l’instant n, avantageusement, mais optionnellement, selon le procédé décrit ci-avant, le procédé se poursuit par l’estimation 504 de leur vitesse moyenne selon ladite direction et à cet instant.
Plus particulièrement, cette estimation débute par une étape 504-2 consistant, pour image binaire redressée IN£, issue du prétraitement de l’image nouvellement acquise INn, à déterminer la distribution DVkopt/n (DV pour « distribution de vitesse ») des pixels de la zone Zkoptn,n (i.e. la zone de l’image IN% associée à la direction Dkoptn) correspondant aux nuages, le long de la longueur L du rectangle Rkopt,n, et donc selon la direction Dk.ovt,n· Ce calcul est par exemple réalisé sur l’image binaire rectangulaire lRkopt,n après rotation de l’angle 90° — akopt en calculant un histogramme comptabilisant, pour chaque ligne de l’image IRkopt/n après rotation, le nombre de pixels de la ligne qui sont égaux à 1, et donc un histogramme selon une direction parallèle à la direction Dk n, tel qu’illustré à la figure 19.
Un flux de distributions FDVkopt§n = [DVkoptil, DVkopti2, -,DVkoptin_ltDVkopt;n^ est ainsi obtenu par l’unité 24, tel qu’illustré à la figure 20 en adoptant la même représentation graphique que précédemment avec l’axe des abscisses correspondant aux instants n d’acquisition. Comme illustré sur cette figure, on constate que ce flux comprend des structures affines qui apparaissent par intermittence. Ces structures correspondent à des passages de nuages. Elles apparaissent d’abord pour des ordonnées hautes (i.e. un banc de nuage apparaît dans le haut de la vignette) puis décroissent en ordonnée au fur et à mesure que ce banc de nuage se rapproche du disque solaire (le centre du disque solaire ayant la coordonnée 0), c’est à dire diminue en ordonnée au fil des images successives dans le flux d’images {INP}. Les inventeurs ont constaté que la pente de ces structures est directement liée à la vitesse moyenne V de déplacement des nuages selon la direction Dkoptin (en pixels par seconde) selon la relation : où β est l’angle que forme la direction des structures affines avec l’axe des abscisses.
Dans une étape 504-4 suivante, le procédé détermine donc la valeur de l’angle β des structures affines en fonction des q dernières distributions calculées {DVkopt:n-p,DVkoptin_p+1,...,DVkopt:n_1,DVkopt^, i.e. sur une fenêtre temporelle [(n — q) x Te; η X Te], où q est entier prédéterminé. La durée de cette fenêtre temporelle est déterminée, par exemple lors d’une campagne de test précédent. Par exemple, une fenêtre temporelle comprise entre 4 secondes et 44 secondes est choisie, ce qui correspond à q= 1, 2,..,11 pour un pas d’acquisition Te de 4 secondes, ce qui a priori est suffisant pour une détermination avec suffisamment de précision de la vitesse moyenne, tout en s’assurant que la vitesse des nuages ne subissent pas de grande modification. D’autres durées sont possibles, par exemple plus importantes. Tout algorithme approprié mis en œuvre par l’unité 24 peut convenir à cet effet (par exemple une régression linéaire).
On note cependant sur la figure 20 que la base des structures (zones proches de la coordonnée 0) peut être plus ou moins large en raison de la proximité avec le soleil. Cet évasement peut perturber la détermination de l’angle par des algorithmes classiques. De manière privilégiée, une approche similaire à celle utilisée pour l’identification de la direction Dkoptn est mise en œuvre. Plus particulièrement, différentes rotations du flux de distribution FDVk n sont testées (correspondant graphiquement à différentes rotations des axes de la figure 20). Normalement, pour une certaine rotation de ces axes, les structures affines deviennent verticales. Dans ce cas, la distribution sur l’axe des abscisses des pixels des nuages présente des pics de valeurs supérieures aux pics des distributions sur l’axe des abscisses des pixels nuages de rotations différentes (les pics sont les plus faibles quand les axes seront tournés d’une valeur conduisant à des structures affines parallèles à l’axe des abscisse). A cet effet, le procédé selon l’invention comprend par exemple : - la production d’une matrice DVn = [DVkoptn_q DVkoptq+1... Wfcopt;„_iWfcopt;„], les distribution DVkoptJl étant codées sous la forme de vecteur colonne de dimension L l’application de différentes rotations d’angles βη prédéterminés et mémorisés dans l’unité 24, et de centre par exemple de coordonnées (n-q, L) à la matrice DVn. Les angles βη sont espacés angulairement entre 0° et 90°, avec un pas angulaire qui règle la précision souhaitée sur la vitesse moyenne recherchée, par exemple un pas égal à 0,5°; le calcul d’un histogramme Histnu pour chaque matrice DVnu obtenue par rotation et le calcul du maximum Mnu de cet histogramme, à savoir le maximum du vecteur [Σ)=ι dvu(l,j) Σ}=ι dvu(2,j) ... Σ)=ι dvu(n-q,j)]; la détermination de l’angle βηορ(,η comme étant l’angle βη correspondant au maximum maxu Mnudes Mn u calculés.
Le procédé calcule alors en 504-6 la vitesse moyenne Vn des nuages, en pixels par seconde, en se fondant sur la relation (2).
La portion supérieure de la figure 2f illustre graphiquement la rotation de la matrice DVn illustrée à la figure 20 pour un angle βη. La figure 22 illustre l’évolution du maximum Mnu en fonction de l’angle βη. On constate que le maximum de la courbe Mnu(/?u) présente un maximum marqué correspondant à l’angle βηορί,η· La figure 23 illustre la détermination manuelle par un opérateur de l’angle βηορΙ,η qui trace les pentes des principales structures de la figure 20. Le rapporteur à droite de la figure 23 montre que l’angle /?optnobtenu par le procédé selon l’invention est cohérent avec une mesure de celui-ci faite directement sur l’image des distributions DVk n.
En outre, l’écart de l’angle /?Uop(.n par rapport à un angle de référence (e.g. 180° si les angles sont mesurés entre 0° et 360° ou 0° si les angles sont mesurés entre -180° et 180°) donne le sens de propagation des nuages selon la direction Dkopt/n. C.5. Prévision de l’ensoleillement
Une fois déterminée la vitesse moyenne des nuages à l’instant n, le procédé se poursuit, en 506, par l’estimation de l’ensoleillement, en particulier l’estimation de la chute d’ensoleillement en raison de nuages passant dans le cône solaire.
En se référant à la figure 24 qui illustre une partie du flux FDVk n, l’axe des ordonnées représente la distance (en pixels) à la base des rectangles de détection (^/copt,n} , c’est-à- dire la distance (en pixels) au centre du disque solaire dans le flux d’images {INn} . Cette distance peut être transformée en temps puisque la vitesse moyenne de déplacement Vn des nuages est connue. Par exemple, en reprenant la vitesse calculée en fonction de la figure 20 (1,48 pixel/seconde), un déplacement de 100 pixels sur l’axe des ordonnées correspond à une durée d’environ 100/1,48 ~ 68 secondes (soit environ 17 images dans le cas particulier d’une période d’acquisition Te de 4 secondes). Autrement dit, une détection de nuage effectuée à l’ordonnée égale à 100 laisse penser que le nuage interceptera le disque solaire 68 secondes plus tard (ou 17 images plus tard). Cela signifie que dans le cas de cette séquence, l’anticipation temporelle maximale est de 500 / 1,48 ~ 340 secondes, soit 5 minutes 40 secondes (500 étant la taille des images redressées et la longueur Z des rectangles de détection).
Toujours en se référant à la figure 24, on peut constater que les nuages en cours de progression vers le disque solaire correspondent aux valeurs significativement non nulles de la distribution DVk n de la dernière image acquise INn. Connaissant donc la position en pixels des nuages à cet instant (i.e. l’ordonnée des valeurs non nulles de la distribution DVkoptlnl ü est ainsi possible de prédire leur moment d’arrivée dans le cône solaire et sur le disque solaire. Le procédé selon l’invention estime donc le passage de nuages devant le soleil en fonction des ordonnées des valeurs non nulles de la dernière distribution DVkopt n et de la vitesse moyenne Vn de déplacement des nuages.
Par exemple, l’estimation de l’ensoleillement consiste à détecter, en 506-2, les valeurs significativement non nulles de la dernière distribution DVk n et à mémoriser leurs ordonnées respectives {Ys} , puis, à chacune de ces ordonnées Ys, à déterminer, en 506-4,
Y l’instant d’arrivée Ts = — des nuages correspondant dans le halo solaire. Les instants
Vn calculés {Ts} correspondent ainsi aux instants pour lesquels il est prédit une chute de Γensoleillement. Les instants {Ts} sont par exemple transmis par l’unité 24 à une unité de contrôle du parc photovoltaïque pour que celle-ci commande en conséquence ledit parc. Une autre utilisation des instants d’arrivée consiste à utiliser ceux-ci pour le diagnostic du bon fonctionnement du parc photovoltaïque. Notamment, un parc est surveillé et des alarmes déclenchées lorsqu’une baisse de production aux origines inconnues est détectée. Connaissant les instants d’arrivées des nuages, cette alarme peut donc être désactivée ou bien classée avec un niveau d’alarme inférieur (e.g. non critique).
Pour illustrer la capacité prédictive du procédé, il a été extrait de l’image de la figure 20 les chroniques du nombre de pixels nuages par ordonnée fixe (c’est à dire correspondant à une échéance temporelle d’arrivée à la position du soleil fixe), à savoir 50, 100, 150, 250, 300 et 350 pixels, tel qu’illustré respectivement aux figures 25A-G. La vitesse Vn étant de 1,48 pixel par seconde, soit 5,67 pixels par pas de temps Te, en décalant temporellement les figures 25A-G de manière correspondante (de 50 / 5,67 = 8,81 pas de temps pour la Figure 25 A, de 100 / 5,67 = 17,6 pas de temps pour la Figure 25B, de 150 / 5,67 = 26,5 pas de temps pour la Figure 25C, de 200 / 5,67 = 35,3 pas de temps pour la Figure 25D, de 250 / 5,67 = 44,0 pas de temps pour la Figure 25E, de 300 / 5,67 = 52,9 pas de temps pour la figure 25F, et 350 / 5,67 = 61,7 pas de temps pour la Figure 25G), les pics des figures ainsi décalées doivent se correspondre et correspondre également à des pics de chute d’ensoleillement. La figure 26 correspond à la superposition des figures décalées. Comme on peut le constater, les pics des différentes figures décalées se superposent, vérifiant en outre que la vitesse des nuages est restée sensiblement constante. La figure 27 illustre un signal enregistré par un pyromètre disposé à proximité du système à des fins de validation du procédé selon l’invention lors de l’enregistrement des images ayant données lieu à la figure 20 ainsi que la superposition des figures décalées. Comme cela est connu en soi, un pyranomètre permet de mesurer la quantité d’énergie solaire incidente sur celui-ci. Comme on peut le constater, les chutes de signal du pyranomètre, qui correspondent à des chutes d’ensoleillement, coïncident avec les pics des figures décalées, prouvant ainsi la capacité du procédé à prévoir l’ensoleillement, et en particulier les chutes d’ensoleillement.
D. CALIBRATION DE LA POSITION REELLE DU DISPOSITIF D’IMAGERIE
Comme décrit précédemment, le positionnement optimal de l’appareil photographique est difficile à obtenir, de sorte qu’il présente une position réelle (0[, 0j,i4zr) différente de la position théorique souhaitée (θ1,θ2,Αζ). Il va à présent être décrit l’étape de calibration 402 mis en œuvre par l’unité 24, ayant pour objectif de déterminer les paramètres de cette correction. Cette calibration peut être réalisée au jour le jour, de manière automatique et/ou sur commande d’un utilisateur, après mouvement du dispositif de prise de vue par exemple (suite à un coup de vent notamment). Cette calibration permet également de détecter si l’appareil a bougé. Notamment, en mettant en œuvre la calibration et en comparant la position réelle obtenue avec une position réelle obtenue au cours d’une calibration antérieure, il est déterminé que l’appareil a bougé si la différence entre les deux positions est supérieure à un seuil. En cas de différence importante, une alarme peut être déclenchée pour obtenir le déplacement d’un opérateur afin de remettre d’aplomb le dispositif de prise de vue.
Dans une première étape 402-2, un flux d’images numériques de calibration {IN£} est acquis par le système 10 alors que le ciel est clair et que le soleil est présent dans la portion de ciel clair photographiée.
Dans une étape suivante 402-4, la position réelle P„ du centre du disque solaire dans chaque image INfi est calculée en fonction de l’image elle-même. Notamment, le calcul par l’unité 24 de la position Pfi du centre du disque solaire met en œuvre un traitement rapide se fondant sur la spécificité de l’image acquise, à savoir celle du ciel, et donc se fondant sur la présence de zones à très forte dominante bleue (le ciel clair) et de zones ayant des tons plus gris (nuage et soleil). Pour chaque image du flux {IN%} , le calcul comprend tout d’abord le calcul de la matrice de luminosité ILcn et de la matrice INn'R~B = |R(IN£) — B(IN£)\. Cette dernière permet de séparer au premier ordre le ciel clair du reste de l’image, les fortes différences correspondant en effet aux parties de ciel clair alors que les faibles différences appartiennent aux nuages et au disque solaire. Un balayage pixel par pixel des matrices Monm 0 IN^R~B et Monm 0 ILcn est ensuite effectué et : si la luminosité du pixel est suffisamment importante (i.e. tZ£(a, b) > su où su est un seuil proche du maximum max{ Monm 0 ILcn) de la matrice Monm 0 ILcn, avec sü compris par exemple entre 80% et 90% de max( Monm 0 ILcn) ou par exemple égal à 250, et que la valeur absolue de r-b du pixel est suffisamment petite, signifiant que le pixel est plus gris que bleu (i.e. ίη£'β-β(α, b) < sr_b où sr_b est un seuil proche de zéro, avec sr_b inférieur ou égal par exemple à 5), le pixel (a, b) est déclaré comme candidat au disque solaire. Les pixels candidats constituent ainsi un sous-ensemble de l’image Monm 0 ILcn; - Ce sous-ensemble est érodé pour éliminer les groupements de pixels candidats de trop petites tailles et pour favoriser l’identification ultérieure de zones distinctes de pixels candidats. Un filtre circulaire d’érosion de diamètre de sept pixels est par exemple mis en œuvre. Il est ainsi obtenu un ensemble de zones de pixels candidats ; le centre de gravité et la distance du centre de gravité à une position prédéterminée a priori (par exemple la position du centre du disque solaire calculée à l’image précédente ou la position théorique Pf du centre du disque solaire calculé par le modèle astronomique optique notamment) sont ensuite calculés pour chaque zone de pixels candidats ; - la zone dont le centre de gravité est le plus proche de la position prédéterminée est retenue comme correspondant au disque solaire ; et le plus grand cercle inscrit dans la zone retenue est alors calculé, et correspond au disque solaire dans l’image. La position finale réelle du centre du soleil Pf retenue pour les calculs ultérieurs mis en œuvre par le procédé est alors le centre de ce cercle.
La figure 28 illustre un exemple de positions réelles Pf déterminée pour le flux de calibration {INf} .
Le procédé se poursuit, en 402-6, pour chaque image du flux {INf} , par le calcul de la position théorique Pf(01,02,Az) du centre du disque solaire calculée par le modèle astronomique optique, cette dernière est calculée à l’aide du modèle astronomique optique connaissant l’horodatage des images et en fixant arbitrairement les valeurs des trois angles (θνθ2,Αζ). En testant un nombre suffisant de combinaisons de valeurs de ces trois angles (θ±, θ2,Αζ), il est alors possible de déterminer la combinaison qui produit la meilleure comparaison entre les coordonnées Pf calculées du centre du disque solaire par cette méthode analytique et les coordonnées Pf du centre du disque solaire par calculées lors de l’étape 402-4. Par exemple, la valeur (θ[,02,Azr) des angles finalement retenue est celle minimisant la somme des différences absolues entre l’ensemble des coordonnées {Pf} et {Pf}. D’une manière générale, toute méthode d’optimisation, e.g. mise en œuvre par ordinateur, pour résoudre un problème d’optimisation selon la relation suivant convient :
où Ω est un ensemble de recherche, par exemple un hypercube centré sur la position souhaitée par l’opérateur lors de l’installation du dispositif d’imagerie (e.g. (0,0,0)) ou M3, Vi = 1,..., n, Pf est un estimateur de la position du disque solaire dans l’image INf se fondant sur une analyse de l’image INf (e.g. la position Pf décrite ci-dessus), et || || est une fonction de coût quantifiant un écart entre deux vecteurs de M3 (e.g. la norme Ll, L2, infinie, etc.).
La figure 29 illustre le résultat de cette calibration en indiquant par une flèche dans l’image l’orientation de l’azimut Nord (également représenté sur la boussole en bas à gauche), ainsi qu’un niveau à bulle symbolisant par un point l’écart à l’horizontalité du dispositif. Sur cette figure, la courbe à l’arrière-plan correspond aux positions du soleil Pn et la courbe au premier plan correspond aux positions du soleil P£ avec les trois angles égaux à (0[, θ^, Azr). On constate que les deux angles définissant les écarts à l’horizontal du fish-eye sont inférieurs à 10°, c’est à dire des angles en effet très petits, compte tenu de la précision découlant de l’expérience (appareil photographique monté sur un pied prenant appui sur revêtement en gravier). Comme on l’observe, les positions théoriques et calculées précédemment empiriquement coïncident très bien. Cette façon de procéder peut être jugée bien plus précise que des tentatives de positionnement à la main, d’un appareil photographique par exemple, sur le terrain et finalement plus simple et plus précise que la méthode classique s’appuyant usuellement sur des images numériques de damiers réguliers pour caler un modèle de distorsion empirique couplé à une projection de perspective.
Il a été décrit un mode de réalisation privilégié dans lequel la direction et la vitesse moyenne des nuages sont calculées pour chaque image acquise, ou de manière équivalente à chaque instant d’acquisition. En variante, un flux d’images numériques est acquis sur une durée prédéterminée jugée appropriée pour le calcul (e.g. d’une durée comprise entre 4 secondes et 44 secondes avec un nombre minimal d’images, compris par exemple entre 2 et 11) et mémorisé dans l’unité de traitement. Une fois la dernière image du flux acquise, l’unité calcule la direction et la vitesse comme décrit précédemment (découpage en zones, calcul de distributions, analyse de la variabilité sur le flux acquis, etc...) en fonction de toutes les distributions calculées sur l’ensemble de ce flux. Le procédé selon l’invention est par exemple mis en œuvre régulièrement (par exemple toutes les 5 minutes ou à chaque fenêtre temporelle acquise) ou sur demande d’un utilisateur. Ceci permet par exemple de libérer les ressources de l’unité de traitement pour d’autres types de tâches.
De même, des traitements supplémentaires peuvent être appliqués. Par exemple, une détection de valeurs aberrantes peut être mise en œuvre (valeurs induites par exemple par le passage d’un oiseau, d’une feuille, etc...) et les images ou les distributions correspondantes écartées. Le calcul est alors par exemple effectué sur une fenêtre temporelle plus importante pour conserver le même nombre d’images ou sur la même fenêtre temporelle mais avec un nombre de distributions moins important.
METHOD FOR ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE
Field of the Invention The invention relates to the field of image processing for determining the position of the sun in an image of the sky. The invention finds particular application in determining the movement of clouds in a sky. The displacement, characterized by the direction and the speed, thus determined can notably be used for the estimation of the sunshine and the electricity production of a photovoltaic park.
State of the art
The irregular and intermittent nature of renewable energies (solar, hydraulic, etc.) implies significant production constraints. This explains why even today conventional electricity production using non-renewable resources (coal, nuclear, etc.), but totally controllable, is used to offset the drop in production by renewable energy in order to guarantee consumers guaranteed access to electricity. In addition, these constraints and their associated costs are all the more important as the fall in production is not anticipated. This is why in the field of photovoltaic production, particular attention has always been paid to the forecast of sunshine, and hence of photovoltaic production that results from it, and this at any time scale (annual, seasonal, daily variability and intraday).
For example, the joint use of meteorological models and satellite images makes it possible to characterize and predict the evolution of sunshine at time horizons ranging from + lh to + 48h. The advection of satellite images makes it possible to predict cloudiness at the shortest time (up to +5h), while statistical adaptation of meteorological models makes it possible to calculate forecasts at longer intervals.
However, at the intraday scale, the energy production on the photovoltaic parks can fall sharply in a few tens of seconds during the passage of clouds above the implantation site. The recovery is just as sudden. This intermittency is problematic when the energy of the park is used in storage / retrieval mode in association with batteries or other means of renewable energy production (hydroelectric power stations). The servo mechanisms on these relatively short reaction times can however be optimized if the energy production by the photovoltaic panels is anticipated from a few minutes to a few hours. In addition, the cloud cover not only influences the amplitude but also the composition of the global radiation. The presence of clouds significantly modifies the relative importance of diffuse and direct radiation. While in clear sky the diffuse radiation represents at best 10 to 15% of the global radiation, the clouds constitute an additional diffusion factor. In high clouds, most of the global radiation is due to diffuse radiation.
However, the techniques based on meteorological model and satellite images do not make it possible to deal with local phenomena or short-time phenomena or to take cloud cover into account. This is why the approach most often adopted to apprehend and anticipate these phenomena in the immediate vicinity of a photovoltaic park is to acquire images of the sky above the park from the ground, to treat them to determine the cloud race in the sky and try to predict on this basis the sunshine above the photovoltaic park.
However, this approach poses problems both in terms of the material used and the image processing performed on the acquired image stream. For example, acquiring the image of the sky in a single take requires the use of an optics having a very large angle of view, for example of the "fish-eye" type which makes it possible to obtain hemispherical images encompassing the entire sky visible from the shooting position. However, the ability to take very wide-angle photos is accompanied by numerous optical defects: significant geometrical distortions, chromatic aberrations as well as other optical defects which later endanger the image processing and / or distort its image processing. result. In addition to the distortion more difficult to correct than for conventional optics, one of the major drawbacks of this type of optics is that the resolution of the image is very poor when approaching the horizon. This is why sophisticated optics have been designed (eg panoramic catadioptric systems based on hemispherical mirrors, convex mirrors with rotating band, etc.). In addition to their complexity, this type of optics requires the use of mobile mechanical parts (eg a tracker to track the path of the sun or a motor to rotate optical elements), which increases the overall cost of device, and poses more problems of reliability and robustness.
Regarding image processing, the most common techniques to track ("tracking") moving objects in a scene are "Particle Image Velocity"("PIV"). This technique uses a search for resemblance between two temporally spaced images of a short time from which a field of local instantaneous displacement vectors (ie: direction of motion and amplitude of velocity) is extracted. While this technique has advantages, including the ability to distinguish uniform areas from highly variable areas, it has many disadvantages. In particular, a significant part of the image is very often not informed about the vectors of displacements and / or is falsely identified as free of object in displacement. This requires a very fine calibration of the treatment and / or a complex and expensive post-analysis in threshold parameters, and therefore particularly complex to calibrate. Moreover a correct identification of the displacements by this technique supposes weak deformations of the objects (clouds) from one image to another, which imposes taking of images at close interval of time. However, clouds are by nature very changeable form, especially because they are deformed, because new clouds appear or disappear, etc ... Also when the clouds approach the solar disk (case particularly interesting since intercepting the solar disk these clouds affect the solar radiation perceived on the ground), their colorimetric properties are often very strongly modified. In particular, they become clearer in the image and thus put the search for similarity in difficulty. In the end, this type of technique is therefore not very robust for determining the stroke of the clouds and requires significant computing resources for the implementation of the treatment.
In particular, the solar disk in an image is a disruptive element that can put in difficulty the subsequent processing to determine the movement of clouds. In addition, knowledge of the position of the solar disk in a stream of images is an important element for predicting sunshine (fall or increase). The precise knowledge of this position requires, however, a precise knowledge of the position of the image acquisition device. This position is, however, entered manually by an operator and corresponds to the desired position for the latter, and not the actual position of the device.
Expose the invention
The object of the present invention is to provide a method for estimating the position of the solar disk in a digital image that allows automatic and accurate calibration of the position of the image acquisition apparatus. To this end, the subject of the invention is a method for estimating the position of the solar disk in a digital image of a sky acquired by means of an image acquisition device comprising an optical device able to focus a scene observed on a focal plane and a digital sensor arranged in the focal plane of the optics, the estimation of the position of the solar disk being carried out according to an astronomical model of the azimuth and the elevation of the sun in the sky having for parameter a predetermined position (θ [, 02, Azr) of the imaging device, defined by the position of the digital sensor with respect to a horizontal plane and by the position of a sensor axis with respect to an azimuth (Azr) .
According to the invention, the method comprises a step of determining the predetermined position (θ [, θ%, Azr) of the imaging device, comprising: - acquisition of a stream of calibrated digital calibration images {IN £} each comprising the solar disk; for each image ΙΝξ of said flux {IN £}, calculating a position P £ of the center of the solar disk in the image IN £ by analyzing the image ΙΝξ; calculating the position of the imaging device (0,> Azr) as a solution of an optimization problem according to the relation:
where Ω is a search set, and || || is a cost function quantifying a difference between two vectors of R3.
According to one embodiment, the images are in color, and in that the calculation of a position P £ of the center of the solar disk in the image IN £ comprises: the calculation of a brightness matrix ILcn of a matrix IN ^ 'R ~ B = | R (INfi) - B (IN £) \ equal to the absolute value of the difference between blue and red channels of the image IN £; and calculating the position P £ in the matrix ILcn.
In particular, the calculation of a position P £ of the center of the solar disk in the image IN £ comprises: - the calculation of pixel areas of the matrix IN ^ 'R ~ Bcandidats at the solar disk, a pixel (a, b ) being considered candidate if its value Z £ (a, è) in the matrix ILcn is greater than a first predetermined threshold and if its value in ^ R ~ B (a, b) in the matrix INn 'is less than a second threshold predetermined; and selecting a candidate zone as belonging to the solar disk according to its proximity to a predetermined position.
According to one embodiment, the method comprises estimating the average direction of displacements of moving objects in a digital image stream {INn} of a scene comprising sky, the digital images of said stream being provided with a reference axis Are ^ the estimate comprising: - the acquisition of the digital image stream {INn}; for each image INn of the flow of digital images {INn} o the calculation of a segmented image IN% identifying the moving objects in the image INn, o the computation of a set of image zones {Zfcn} respectively associated with analysis directions {Dk} forming different angles {ak} predetermined with the reference axis Aref: o for each zone Zkn of said set of areas {Zkn}, calculating a distribution DDkn of the values of pixels of said segmented image IN £ included in said zone Zkn, the distribution DDkn being calculated along a direction normal to the analysis direction Dk of said zone Zkn; for each analysis direction Dk, the computation of a variability criterion CVkn of distributions {DDkn_p, DDkn_p + 1 ,, DD ^ n ^, DDk> n} calculated for zones Zkn associated with said direction Dk; and determining the average direction Dkoptn of moving objects in the digital image stream {INn} as being equal to the analysis direction for which the variability criterion CVk n is optimum.
In other words, the inventors have discovered that the oriented zone which has the minimum temporal variation of its distribution is that whose direction corresponds to the average direction of movement of the objects in a scene. This process is well adapted to the race of clouds in the sky. It is also robust because it naturally filters the phenomena of high dynamics / short duration which usually disturb the treatments of the state of the art. In addition, except possibly the segmentation of the images, the calculations implemented are basic (calculation of distributions and calculations of distribution variations), and therefore do not require large computing resources. In addition, if the segmentation may require resources or a calibration, the method according to the invention is however very robust to it so that a low precision segmentation of the images is sufficient. Finally, it has been found that the method is also robust to the optical process of image formation and returns very satisfactory results, even on images made from a fisheye, including when applied to them a rectification without treatment of aberrations related to this type of optics.
According to one embodiment, the acquisition of the digital image stream {/ n} comprises: the acquisition of the scene by a shooting device having a field angle greater than or equal to 100 ° according to the diagonal of a 35mm format or equivalent 35mm format sensor, in particular a 180 ° angle of view; and prior to the calculation of the DDkn distributions, the straightening of the images the digital images {INn} The invention thus advantageously exploits the optics with a very angle of field, for example the fisheyes which make it possible in a single shot to acquire the all of the visible scene, and therefore the visible portion of the sky.
According to one embodiment, the scene is acquired using an imaging device comprising an optical system capable of forming an image of the scene on a focal plane and a digital plane sensor disposed on the focal plane of the optical system. for producing the digital image stream {INn}, and the method comprises computing the set of image areas {Zkn} for each INn image of the digital image stream {INn}, said calculation comprising: computing a set of rectangles {Æfcn} each having a small side centered on the same predetermined point of the focal plane, and having their long sides respectively parallel to the analysis directions {Dk}; and calculating the areas {Zfcn} as being respectively equal to the intersection of the rectangles {Æfcn} with the image INn.
In particular, the scene comprises sky, and the predetermined point of the focal plane is the position of the center of the solar disk on the focal plane of the optical system.
In other words, the method makes it possible to directly calculate the average displacement in a reference frame of interest, for example centered on the image of the sun. It is thus determined directly clouds that move to the sun and are therefore likely to hide the latter later.
Advantageously, the small side of each rectangle Rkn of the set of rectangles {Æ / c, n} is at least equal to the diameter of the halo of the solar disk. The images are thus divided into zones which directly take into account a result sought after, notably the occultation of the solar disk and its halo by clouds and thus the variation of the sunlight.
According to an advantageous variant, the position of the center of the solar disk is calculated according to an astronomical model of the azimuth and the elevation of the sun in the sky defining the position of said center on the focal plane. In other words, it is possible to determine the movement of the clouds towards the sun, even when the latter is absent from the image, for example when it is hidden by a cloud, by a tree, or even completely by outside the image.
According to a particular embodiment, the variability criterion CVkn for the direction Dk is calculated as a function of the distributions {DDkn_p, DDkn_p + 1, DDkn] associated with the last p> 2 digital images acquired.
According to one embodiment, the variability criterion CVkn for the direction Dk is calculated according to the relation:
an expression in which j refers to the position of a digital image in the digital image stream {INn}, n and np respectively refer to the nth image and the (np) th image in the digital image stream {INn}, with p strictly less than n, k referencing the kieme direction Dk, and || L || is a norm of the vector V of the space M1, / being the dimension of the distributions DDkn.
The inventors have indeed noted that the CVkn criterion has an overall minimum as a function of the index k and that this global minimum corresponds to the mean direction of displacement. The criterion can be more or less refined to emphasize this overall minimum, but even the simple sum of the two-by-two distribution difference values makes it possible to detect this minimum with precision.
According to one embodiment, the method comprises calculating an exclusion mask of the immobile elements of the scene, and prior to calculating the DDkn distributions, said mask is applied to the images of the digital image stream {INn} and / or segmented images {/ N%}. In other words, excluding objects in the scene that do not contain any significant information (trees, reliefs, buildings, etc.), the treatment focuses only on the sky and thus sees its precision reinforced.
According to one embodiment, the method further comprises: calculating an average speed Vn of moving objects in the digital image stream {INn} according to said calculated mean direction of displacement; and predicting the position of an object in motion according to said average velocity Vn and said average direction Dkoptn.
In particular, moving objects are clouds, and the method includes cloud-based sunshine forecasting in the digital image stream {INn} having said velocity and said mean direction. In particular, the method includes predicting sunshine over a photovoltaic park, in particular the fall or increase of sunshine.
In other words, the method advantageously exploits the fact that the average direction of displacement is already known for calculating the average speed in this direction. Having an important information, the calculation of the speed is thus facilitated, any type of calculation can be used for this purpose. Knowing the direction and the average speed, for example clouds in the sky, it becomes possible to make intra-day forecasts of the sunshine, and therefore of the electricity production by a photovoltaic park.
According to one embodiment, the calculation of said average speed V n comprises: for each segmented image IN%, the calculation of a distribution DVk n of the values of the pixels of the segmented image IN £ belonging to the zone Zk n associated with the determined mean direction DkoptTl, the distribution DVk n being calculated along the determined mean direction Dkopt Tl; calculating the inclination of at least one affine structure included in a matrix DVn = [DVkoptn_qDVkoptn_q + 1 ... DVk () ptn_iDVkopt1 ^ resulting from the concatenation of the DVk distributions n of the images between the (nq) th image of the stream up to the nieme of the digital image stream {INn}, where q is an integer greater than or equal to 1; and calculating the average speed Vn as a function of the calculated inclination.
In other words, the inventors have found that naturally occurring affine structures in the matrix DVn = [DVkoptin_qDVkoptin-q + 1 -DVkoptin-if> Vfcopt, n] · These affine and parallel structures correspond to the average displacement of the objects in the image flow and their slope corresponds to the average speed of the objects according to the average direction of movement. The calculation of the average speed is thus robust and is limited to the determination of a slope.
According to one embodiment, the calculation of the inclination of the at least one affine structure included in a matrix DVn comprises: the application of different rotations of predetermined angles βη to the matrix DVn; calculating a histogram Histnu for each matrix DVnu obtained by rotation and calculating the maximum Mn u of this histogram; the inclination determination being equal to the angle βηορ (, η corresponding to the maximum of the maxima Mn u calculated.
This calculation makes it possible to accurately determine the inclination using simple algebraic calculations, and therefore little consumer of computing resources. The invention also relates to a system for estimating the position of the solar disk in a digital image: a camera for acquiring the digital image; an information processing unit configured to estimate the position of the solar disk in the image as a function of an astronomical model of the azimuth and of the elevation of the sun in the sky having as parameter a predetermined position (0 , 0J> Azr) of the imaging device, defined by the position (0 [, 0 £) of the digital sensor with respect to a horizontal plane and by the position of an axis of the sensor with respect to an azimuth (Azr); and said unit being further configured to determine the predetermined position (0 [, 6r2, Azr) of the imaging device: for each image ΙΝξ of a timestamp digital calibration image stream {IN £} each comprising the solar disk {IN £}, by calculating a position P £ of the center of the solar disk in the image IN £ by analyzing the image IN £; computing the position of the imaging device (0j \ 0j, Azr) as a solution of an optimization problem according to the relation:
where Ω is a search set, and || || is a cost function quantifying a difference between two vectors of M3.
In particular, the system is configured to implement a method of the aforementioned type. The invention also relates to a computer program product recorded on a computer-compatible computer medium comprising instructions for the execution of a method for estimating the position of the solar disk in a digital image of sky acquired by means of a computer. an image acquisition device comprising an optical device capable of focusing an observed scene on a focal plane and a digital sensor disposed in the focal plane of the optics, the estimation of the position of the solar disk being carried out according to an astronomical model of the azimuth and the elevation of the sun in the sky having for parameter a predetermined position (θ [, θζ, Αζ1 ") of the imaging device, defined by the position (θ [, θ £ ) of the digital sensor with respect to a horizontal plane and by the position of a sensor axis with respect to an azimuth (Azr), the method comprising a step of determining the predetermined position (β ^, θ ^, Azr) of d imaging device, comprising: acquiring a stream of calibrated digital calibration images {ΙΝξ} each comprising the solar disk; for each image IN £ of said flux {IN £}, calculating a position P £ of the center of the solar disk in the image IN £ by analyzing the image ΙΝξ; calculating the position of the imaging device (0,> Azr) as a solution of an optimization problem according to the relation:
where Ω is a search set, and || || is a cost function quantifying a difference between two vectors of M3.
In particular, the support comprises instructions for the implementation of a method of the aforementioned type.
BRIEF DESCRIPTION OF THE FIGURES The invention will be better understood on reading the following description, given solely by way of example, and with reference to the accompanying drawings, in which: FIG. 1 is a diagrammatic view of a sunlight prediction system according to the invention; FIGS. 2A-2D are flowcharts of a prediction method of prediction of solar photovoltaic production implemented, for example implemented by the system of FIG. 1; Fig. 3 is an image, transformed into black and white, of an image acquired by the system of Fig. 1; FIG. 4 is an image derived from a segmentation, eg of Otsu type, of FIG. 3 for the determination of non-meteorological objects; Figure 5 is a profile illustrating the location of a portion of sky in the image of Figure 3, from the segmentation of Figure 4; Fig. 6 is a digital image of an exclusion mask obtained from the profile of Fig. 5; FIG. 7 is a digital image of an exclusion mask obtained by applying a circular expansion mask in the image of FIG. 6; FIG. 8 is a digital image illustrating the application of the mask of FIG. 7 to the digital image of FIG. 3; FIG. 9 is a digital image illustrating the application of a mask further comprising the solar disk to the digital image of FIG. 3; FIG. 10 is a digital image illustrating a segmentation for the detection of clouds in the image of FIG. 9; FIG. 11 is an image illustrating a detection rectangle according to the invention; FIGS. 12A-C are figures illustrating different type of detection zones according to the invention; FIGS. 13A-E are images illustrating flow of displacement distributions for different angles of detection; FIG. 14 is a schematic view graphically illustrating the calculation of a distribution flow; FIGS. 15 to 17 are plots illustrating different criteria of flux variability as a function of the detection angle; FIG. 18 illustrates a comparison of the automatic detection of the average displacement direction according to the invention and a manual detection of this direction; Figure 19 is a schematic view graphically illustrating the calculation of a speed distribution according to the invention; Fig. 20 is an image graphically illustrating affine structures of a stream of velocity distributions; Fig. 21 is an image illustrating a rotation of the structures of the figures of Fig. 20; FIG. 22 is a plot of a criterion of a speed distribution as a function of an angle of rotation of this distribution; FIG. 23 illustrates the comparison between the automatic determination of the angle of rotation corresponding to the speed of movement of the clouds and the manual determination of this angle; Fig. 24 is an image illustrating a portion of the speed distribution flow; FIGS. 25A-G are plots illustrating temporal arrival times of cloud on the solar disk for different distances in pixel from the sun, FIGS. 26 and 27 are plots superimposing the estimated times of arrival of the clouds on the sun, FIG. 27 further illustrating measurements of illumination by a pyranometer; and Figures 28 and 29 are images illustrating the method of calibrating the position of the camera according to the invention.
Detailed description of the invention A. Ratings
In the rest of the document, the classical mathematical conventions of linear algebra are applied for real matrices of M / 1xS, where A and B are positive integers, with in particular for matrices M and N of RAxB: M - (ma , ù) ae [1A] jùe [1B] \ M \ - (\ ma, b \) ae [lA], be [l, B]
The notation Μ ® N relates to the application of an exclusion mask Ma the matrix N. More particularly, the application of a treatment on the object M ® N means that only the pixels of the corresponding matrix N M mask pixels having a first value are taken into account for processing. For example my b = 1 if the pixel is taken into account or 0 if it is excluded from the treatment. Similarly, when reference is made to the object Μ ® N alone (the "matrix M ® JV"), reference is made to the set of pixels of N which are not excluded by the mask M. to take or not account a pixel of the matrix N according to the mask M can be achieved in many ways known in itself, the concept of mask being conventional in image processing. By convention, in the application, the "matrix" M ® N is graphically represented by the image of N modified by setting its pixels to black when they correspond to the pixels of the mask M of value 1. B. MATERIAL ·
Referring to FIG. 1, an imaging and image processing device 10 according to the invention, directly implanted on a photovoltaic park, comprises a camera or a camera 12, for example of the reflex type, equipped with a CMOS or CCD digital sensor 14 for the acquisition of an image in the visible radiation, for example in the 35mm format, and a lens 16 with a very wide field angle of greater than 100 ° on the diagonal of the sensor 14, advantageously 180 ° (eg fisheye). The objective thus produces on the sensor 14 a circular optical image of the space in front of the objective 16 and therefore an image comprising all of the sky visible from the position of the device 10. The sensor 14 is arranged on the focal plane 18 of the objective 16 and centered on the optical axis 20 of the latter. The apparatus 12 further comprises the circuits 22 necessary for controlling the acquisition of the images and the pre-processing of the data from the sensor (eg demosaicing) to produce an image from the signals generated from the sensitive photosites of the sensor 14, usually a "DSP" processor (for "digital signal processor"), as well as a clock for timestamping acquired images Acquisition of images is performed with an acquisition period Te, the device comprising an internal intervalometer for this purpose or the device having an external intervalometer in a manner known per se.For example, the apparatus 12 is a Canon 6D reference camera equipped with a Canon EF 8-15mm f / 4L lens set to the 8mm focal length, piloted by a DigiSnap 2700 external reference intervalometer from Harbonics.
The device 10 furthermore comprises an information processing unit 24 connected to the apparatus 12, by a cable or by a wireless communication, for receiving the images acquired by the apparatus 12, for example in raw format (after the demosaicing but before any computer processing, corresponding to the format "RAW" at Canon) or in a reading format (eg jpeg), and for applying a processing implemented by computer according to the invention. The unit 24 is for example a personal computer, a tablet, a smartphone, or more generally any microprocessor-based system s), particularly of the DSP ("digital signal processor") type, based on FPGA-type circuits, base of circuits mixing these types of technology, etc., able to implement the invention. The unit 24 is notably provided with all the memories (RAM, ROM, cache, mass memory, etc.) for storing the instructions implementing the method according to the invention, parameters and settings that are useful for this purpose. implementation and memorization of the results of intermediate and final calculations.
Optionally, the camera 12 is advantageously housed in a protective case 26 provided with a transparent window 28 in front of the lens 16 to protect it from the weather and thus allow the device to remain permanently on the chosen site. In a variant, the unit 24 is also housed in the protection box 26 and communicates, via for example a wireless network (telephone, satellite, etc.), the result of the process that it implements to a remote site where the production of the photovoltaic park is managed.
In operation, a stream of timestamped digital images {INn} is thus produced by the apparatus 12, a digital image INn being acquired for example every Te seconds, the integer n representing the nth acquisition time nxTe for a moment initial n = 1, and therefore the // 'th picture in the image stream {INn}. In a variant, the acquisition is performed in a non-regular manner with a non-constant acquisition step in time.
Each INn image comprises A lines of pixels per B columns of pixels and comprises three channels of a color space, for example the red (r), green (v) and blue (b) channels respectively coded by three matrices R, F, RAxB B. The numerical values composing the matrices are coded on a predetermined number of bits, for example 8 bits, corresponding to 256 levels of values (from 0 to 255). Subsequently, the channels r, v, b of any color image "IN" are respectively denoted R (IN), V (IN) and B (IN).
In what follows, the color space considered will be the RGB space, it being understood that the invention applies to other color spaces, for example the HSV space (hue, saturation, light intensity). Similarly, a number of coding bits equal to 8, corresponding to values between 0 and 255, is described. Naturally, the invention applies to any number of bits. C. Computer processing
It will now be described, in connection with the flow charts of FIGS. 2A-2E, a method according to an embodiment of forecasting the electricity production of a photovoltaic park, and more particularly the fall of sunshine on the park. photovoltaic causing a fall in production, implemented by the device of Figure 1, this method based on the estimate according to the invention of a mean direction of movement and an average speed of the clouds in the sky. Cl Calibration
The method starts with a phase 40 of positioning and calibration of the imaging device 10.
A first step 400 consists in loading into the memory of the unit 24 all the instructions and parameters used by the method, in particular: an astronomical model of the position of the sun in the sky, for example determining the azimuth and the elevation of the sun in the sky according to the longitude, latitude, time and day of the year. This model is for example the "solar position calculator" of the NOAA ("National Oceanic and Atmospheric Administration"); a model of the image formation on the sensor 14 by the objective 16. For fisheye, it is for example the model described in the document "Geometry modeling and calibration of fisheye lens camera System", Schwalbe E. In Proc. 2nd Panaramic Photogrammetry Worshop, Int. Archives of Photogrammetry and Remote Sensing (Vol 36, No. Part 5, W8), February 2005. This model makes it possible, in particular, to model the image of an object photographed by the device 12 in the focal plane of the objective 16 , and thus on the sensor 14.
By knowing the position θ2, Αζ) of the sensor 14 in a terrestrial reference, and therefore the position of the sensor relative to the sun, the date and the time, the combination of the astronomical model and the model of the image formation makes it possible to calculate the image of the center of the sun on the focal plane of the lens, and therefore the center of the solar disk in the images produced by the camera 12. This combination is hereinafter referred to as the "optical astronomical model".
In 402, an operator positions the apparatus 12 at a selected site with the sky objective. The position of the apparatus is thus defined by three parameters (01 # θ2, Αζ), the pair of angles (θ ±, θ2) defining the position of the sensor 14 with respect to a horizontal plane, and the parameter Az representing the angle that forms the axis passing through the centers of the top and bottom edges of the sensor with an azimuth, for example the north azimuth. The operator seeks in particular to position the sensor 14 on a horizontal plane with the north azimuth corresponding to the middle of the upper edge of the images produced by the apparatus 12. The desired, ie ideal, position of the apparatus is then equal to ( θνθ2, Αζ) = (0,0,0). The geographical coordinates, eg longitude and latitude, of the device are also stored in the unit 24.
Obviously, it is desired to position the device in an open area so as to maximize the sky portion in the images. However, the device can be implanted in a very constrained zone in terms of irrelevant elements present in the field of the device (or "non-meteorological object", or "onm", for example trees, constructions, pestles , etc ...), the invention allows indeed a robust prediction including small portions of the sky and including the sun hidden.
Once the device 10 is positioned, a calibration 404 is implemented. In particular, the perfect positioning of the apparatus is difficult to obtain, so that it has a real position (θ [, 02, Azr) different from the theoretical position (θ1, θ2, Αζ). The position of the device being used during processing steps according to the invention, a correction of the difference between the real and theoretical positions is implemented. The purpose of the calibration step 402 is to determine the parameters of this correction and is described in more detail below. The corrected position (θ [, 02, Azr) is memorized by the unit 24 and is that used by the optical astronomical model to calculate the center of the solar disk in the images produced by the apparatus 12.
The initial phase 40 is continued by calculating, at 406, an exclusion mask of the non-meteorological objects of the scene present in the acquired image stream. This mask makes it possible to define the pixels of the images on which the treatment will be applied for the forecast of the sunshine. Being in a phase of initialization and calibration, it is possible to implement any known method to determine this type of mask, including the most sophisticated methods and consuming time and material resources, or even with the operator assistance who designates by hand all or part of these objects.
Advantageously, the determination of the exclusion mask is based on the specificities of the images acquired namely that of the sky, which allows an automatic and rapid determination of the mask. Notably, during the day, the "clear sky" (portion of cloudless sky and without the solar disk, ie the portions of blue sky), the clouds, the solar disk and the non-meteorological objects have different luminosities, as illustrated for example in Figure 3.
The determination of the exclusion mask thus comprises: a) the acquisition by the device 12 of a digital image INm during the day, for example in sunny or slightly cloudy weather with the sun high in the sky, eg at its zenith; b) calculating a brightness matrix ILm of the image INm, in a manner known per se, for example according to the formula ILm = -X (R + V + B), thus obtaining a gray-scale image as shown in Figure 3; c) determining a first mask Mx G RAxB corresponding to the optical image of the scene on the sensor, namely m1 (a, b) = 1 if the pixel (a, b) of INm belongs to the image of the scene on the sensor 14 (ie to the fisheye disk on the image) and 0 otherwise; d) the application of a brightness histogram segmentation algorithm, as for example described in the document "A threshold selection method form gray-level histograms", by Otsu N, Automatica, 11 (285-296), 23-27, 1975. This type of algorithm consists in dividing the histogram of the matrix Mt 0 ILm into a plurality of homogeneous parts from the point of view of their values, for example four parts by calculating three thresholds SI, S2, S3 such that: (51.52.53) = argmin (ff (0, sl) + a (sl + 1, s2) + a (s2 + 1, s3) + a (s3 + 1.255)) (sl, s2, S3) e [0,255] 3 where σ (η1, η2) is the standard deviation of the values of the abscissa histogram between ni and n2. By adopting 5 gray levels NG1 = 0 <NG2 <NG3 <NG4 NG5 = 255, the resulting image M1 0 of this segmentation (il ^ tu {a, b) = N G1 if ilm (a, b) <51, you have (a, b) = NG2 if 51 <ilm <52, etc.) is illustrated in Figure 4. This type of problem is simple resolution using a suitable solver, for example implemented by the unit 24; e) For each of the angles of a trigonometric circle centered on the center of the image (between 0 ° and 360 ° with a predetermined angular pitch, for example of 1 °), the pixel of the matrix ILm is determined. closer to the center of the circle whose brightness iZm (a, b) is below the threshold SI which defines the darkest brightness class resulting from the segmentation. It is thus obtained a profile, for example stored in an array of 360lp values of angles, of the distance to the center of the image which is the smallest distance identified of a dark pixel, the resulting profile of the image Mt 0 ILm being illustrated in Figure 5; f) the determination of a second mask M2 e Ey4xBwith m2 (a, b) equal to 0 for all the pixels beyond and on the profile and equal to 1 otherwise, the mask M2 being illustrated in FIG.
The mask M2 is in a first variant that used during the subsequent treatment of sunshine forecasting. However, in Figure 6, we see that the limits of the masked area (in black) are very close to the tree crowns and the contours of the branches. One suspects that in case of wind, this mask would not exclude the branches that would move and appear then in the visible part (white) of the mask. In a second variant, an additional correction of the mask M2 is applied to take account of this phenomenon. In particular, an expansion algorithm is implemented on the mask M2 with a radius sufficient to encompass at least the foliage, for example by means of a circular expansion mask having a radius of between 10 and 50 pixels. More particularly, for each of the pixels previously defined as masked of the mask M2 (m2 (a, b) = 0), a circular mask is centered on this pixel and all other pixels of M2 contained in the circular mask are set to 0. The M3 mask resulting from the expansion of the M2 mask using a circular expansion mask of radius equal for example to 10 pixels is illustrated in Figure 7. Referring to the representation of the matrix M3 0 ILm of the figure 8, it can thus be seen that the exclusion mask makes it possible to identify the useful parts of the image (here the sky, the cloud and the sun) on which the subsequent treatments implemented in the context of the prediction of the fall of sunshine on the photovoltaic park.
The exclusion mask corresponding to the non-meteorological objects Monm = M2 or Monm = M3 is then stored in a memory of the unit 24.
Once the initial phase of positioning and calibration is complete, the device is put into service. The apparatus 12 produces the digital image stream {INn} which is processed at 50 by the information unit 24 for estimating sunshine. Since the method is intended to determine the intraday variations in sunlight, the complete treatment is advantageously implemented only when there is clear sky. For example, the unit 24 determines the average brightness of each image INn and the unit 24 triggers the processing if the average brightness (eg of the image Monm 0 ILm) exceeds a predetermined threshold, meaning that the sky is not covered. Other tests, or no test, trigger can of course be implemented. The computer processing step 50 comprises, for each new digital image acquired INn, a step 500 of preprocessing the image, followed by a step 502 of estimating the mean direction of the clouds at the acquisition instant n of the image, followed by a step 504 for estimating the average speed of the clouds at instant n, followed by a step 506 for predicting the fall of the sunlight above the system 12, and therefore of the park where this the last one is implanted, according to the direction and the average speed computed at the time of acquisition n. C.2. Pretreatment of the image lNn
In order to facilitate the subsequent segmentation between clear sky and clouds, the process starts, in 500-2, with the additional exclusion of the solar disk if it is present at least partially in the new acquired INn image, which leaves no subsist that the portions of the acquired image INn corresponding to the clear sky and the clouds. More particularly, the unit 24 calculates the theoretical position Pn of the center of the solar disk in the focal plane and therefore, if appropriate, in the image INn as a function of the optical astronomical model, then determines a disk of radius rs having as its center the position Pn. The radius rs, parameter recorded in the unit 24, is a value determined for example during a previous test campaign and is greater than or equal to the largest radius of solar disk observed during the year, or a defined period of the year, over the whole year or over several years. Other methods of calculating said radius are obviously possible. The ray rs thus allows to exclude in a certain way the image of the sun in the digital image INn without excluding a too large area. This value, which depends on the resolution of the image, is for example equal to 80 pixels.
An exclusion mask M4n of M / 1xB is then calculated by the unit 24. This mask is equal to the mask Monm when the solar disk is not included in the image. In the opposite case, the new exclusion mask M4n corresponds to the exclusion of the non-meteorological objects and the solar disk pixels included in the INn image, namely m4 "(a, b) = 1 if m3 (a, b) = 1 or if (a, b) belongs to the solar disk, and m4n (a, b) = 0 otherwise. The image M4n 0 ILn thus corresponds only to the portions of clear sky and clouds of the photographed scene and is illustrated in FIG. 9.
A segmentation of the image M4n 0 INn is then implemented, in 500-4, to identify in it the portions of clear sky and clouds.
Advantageously, the matrix M4 n 0 INB ~ B is used for this segmentation, where IN * ~ b = | R (/ Nn) - B (jNn) |. As described below in relation to the calibration of the actual position of the apparatus 12, the matrix INB ~ B, resulting from the difference of the red and blue channels, shows in a contrasting manner the portions of clear sky and the clouds. A simple threshold on the central value of the histogram of this matrix can thus suffice to segment the image. An algorithm of the Otsu type as described above, with a single threshold, is for example alternately used to distinguish in the histogram the class of the clear sky and the class of clouds. This allows a more precise segmentation, especially concerning the portions of clouds. Following this segmentation, a binary matrix M4 n 0 IN £ cn is thus calculated ("CCN" for "clear sky / cloud") with for example in ^ CN (a, b) = 1 if the pixel (a, b) belongs to a cloud and 0 if the pixel (a, b) belongs to a portion of clear sky. This matrix is illustrated in Figure 10.
The pretreatment 500 is continued at 500-6 by the correction of the geometric distortion induced by the objective 16, in particular the correction of the masks M4n and the binary image INfiCL. This type of treatment, also known as "straightening", aims for a straight line in the photographed scene to be also a straight line in the straightened image.
In a preferred manner, the recovery implemented is based on the fact that the invention is insensitive to the effects of perspectives in the image, and therefore advantageously uses the simplified recovery of the projection of a 3D scene on a 2D plane by a fisheye, for example that described in the Schwalbe document E. cited above This rectification, which only corrects the image deformations related to the projection through the lens, is very simple and quick to calculate.
In the following, the clear sky portions being no longer used as such for determining the direction of clouds, the method continues for example by generating a rectified binary image, denoted IN £, identifying the clouds INn image, ie in £ (a, b) = 1 if the pixel (a, b) belongs to a cloud and 0 otherwise. This image therefore includes the mask M4 n.
Finally, the unit 24 calculates at 500-8 the position Sn = (xn, yn) of the center of the solar disk in the rectified image IN £. The pretreatment implemented by the unit 24 thus calculates for each digital image INn acquired by the device 12, a binary image rectified IN%, and the position Sn of the center of the solar disk in the rectified image.
Note that the dimensions of the rectified images are not necessarily equal to those of the images from the sensor. For the sake of brevity, however, this assumption will be made later, since the matrices and calculations can of course be defined for different dimensions. The same hypothesis is used for the rotation of a matrix.
The method continues by estimating the average stroke of the clouds at time n, and more particularly by estimating 502 of their mean direction, followed by the estimate 504 of their average speed in the mean direction. C.3. Estimation of the average direction of the clouds The estimation of the mean direction of the clouds consists in a first step in cutting the new image IN £ resulting from the preprocessing of the image INn into a plurality of zones, each associated with a particular search direction Dk predetermined and constant in the stream of acquired images {INn}. Then, in a second step, the estimation determines the evolution of the clouds as a function of time in each zone. The direction Dk of the minimum cloud variation zone then corresponds to the average direction of the clouds. Advantageously, the estimate according to the invention is configured to take account of the desired sunlight forecast.
Referring to FIG. 11 illustrating a rectangle on a rectified M4n®INn image, step 502 for estimating the mean direction of clouds begins, in 502-2, with the definition of a set of detection rectangles { /? fcn} of length and width in identical pixels but of different orientations. More particularly, the detection rectangles Rkn have: a width in pixels equal to at least the diameter of the solar halo, denoted D%. The solar halo, also called "solar cone", corresponds to a zone of luminosity very close to that of the solar disk and surrounding it, whose occultation by clouds thus also induces a fall of the sunshine. The width / is for example determined during a previous test campaign, stored in the unit 24 and is greater than or equal to the maximum diameter of the solar halo observed during this campaign. For example, for an IN £ image of 500 pixels by 500 pixels, the width / is equal to 160 pixels; - a small side centered on the center Sn of the solar disk. The rectangle thus makes it possible to take into account all the clouds that may later intercept the solar halo D%, and thus vary the sunshine; a length in pixels L, for example equal to the largest dimension of the images IN £ (eg 500 pixels in the examples illustrated in the figures). As detailed below, the length of the rectangles defines the temporal window for predicting the variations of sunshine. In particular, the greater the length L is important and the more distant clouds of the sun are taken into account for the forecast of sunshine, in particular the fall of sunshine; a predefined direction Dk, parallel to the length L, indicated by an angle ak with respect to a predetermined axis of the image, for example a horizontal axis of the images (eg their lower edge). The directions Dk are for example evenly spaced angularly between 0 ° and 360 °, with an angular pitch Θ that sets the desired accuracy on the desired average direction. For example, step Θ is chosen equal to 1 °, resulting in K = 360 fixed Dk directions studied.
As illustrated in FIGS. 12A-C, depending on the position of the solar disk in the rectified image, a rectangle may be included entirely in the image (FIG. 12A), or partially in the image if the solar disk is in the image. image and close to an edge thereof (Figure 12B) or if the solar disk is out of the image and close to an edge of the image (Figure 12C).
The method consists, in a next step 502-4, of determining, for each zone Zfc n of the rectified binary image IN £, included in a rectangle Rkn, the distribution DDkn of the pixels (DD for "displacement distribution") corresponding to the clouds along the width / rectangle Rk n, and therefore in a direction normal to the direction Dk associated with this rectangle.
For example, for each detection rectangle Rk n, the method: produces a rectangular binary image lRkiU corresponding to the content of the rectangle Rkn (ie the zone Zfcn, where appropriate supplemented by a portion of pixels at 0 corresponding to the area of the rectangle Rkn not included in the IN £ image, as shown in Figures 12B and 12C); applies 90 ° - ak angle rotation to obtain a rectangular image IRkn of lx L pixels, the small side of the rectangle thus being horizontal; calculates the DDkn distribution. For example, the distribution DDkn is a histogram counting for each column of the image lRkn the number of pixels of the column which are equal to 1, ie corresponding to clouds.
For each digital image INn acquired by the device 10, the unit 24 thus calculates a set of distributions {DDkn} respectively associated with fixed directions {Dk}. The unit 24 therefore produces, for the image stream {INn}, K stream of distributions, for example noted for the eighth stream FDDkn = {DDk1, DDk2, -, DDkn_1, DDkn}. Examples of FDDkn distribution flow are illustrated graphically in FIGS. 13A-13J which correspond respectively to the angles ak equal to 0 °, 10 °, 20 °, 30 °, 40 °, 50 °, 60 °, 70 °, 80 ° and 90 °. Referring to FIG. 14, the graphical representation consists in transforming the scale of the values of the distributions (from 0 to I concerning histograms) into a scale of gray, which gives a graphical representation of a distribution in the form of a label, then juxtapose the labels according to the increasing values of n.
According to the invention, the analysis of the variability of the FDDkn distribution flows makes it possible to determine the mean direction of the clouds at time n. For this purpose, in a next step 502-6 of the method, the variability at instant n is calculated as a function of the last p calculated distributions {DDkn_p, DDfcn_p + 1, ..., DDkn_v DDkn], ie over a time window. [(n - p) xTe; nx Te \, where p is a predetermined integer. The duration of the window is for example between 4 seconds and 44 seconds, which corresponds to p = 1, 2, .., or 11 for an acquisition step Te of 4 seconds. The variability is for example calculated by any appropriate criterion CVkn, in particular a criterion of resemblance or dissimilarity between the distributions comprised between the moment of acquisition n and the previous moment of acquisition np.
Advantageously, the criterion CVkn is equal to the sum of the resemblances between two distributions spaced by m> 1 time step, namely a CVkn dissimilarity criterion for m = 1 according to the relation:
(1)
The inventors have discovered that such a criterion comprises two optimums (two maxima if similarity criterion or two minima if dissimilarity criterion) as a function of the index k. These angles correspond respectively to the direction of the angle of the clouds approaching the solar disk and the angle of the cloud pixels moving away from the solar disk (ie corresponding to an angle ak and an angle ak + 180 °). In other words, if the clouds come from the northwest, the resemblance of the distributions is maximum for a detection angle akn of 135 ° but also for the angle 135 ° + 180 ° = 215 ° (south-east direction) since, according to this direction, the same clouds which, once passed under the solar disk, move away with the same direction.
The method thus determines the index k, corresponding to the angle ak, which maximizes (resp., Minimizes) the criterion CVkn between 0 ° and 180 °, as well as the index k \ corresponding to the angle ak + 180 °, which maximizes (resp. minimizes) the criterion CVkn between 180 ° and 360 °. For example, for an angle step Θ of 1 °, the method determines the index k = kopt between 0 and 180 which maximizes (respectively minimizes) the criterion CVkn = CVkn + CVk + 18Q n. The direction Dkoptn corresponding to the index kopt determined by the unit 24 for the moment n, is thus the mean direction of movement of the clouds at time n.
Figure 15 illustrates CVkn criteria based on relation (1) for different temporal distances between distributions. The circles indicate the value of the angle (Xkopt For which the criterion CVkn is minimal, and therefore the direction Dkoptn of propagation of the clouds at the instant n finally identified by the method according to the invention. sums is weakly sensitive to the temporal gap between two distributions.In particular, the cloud propagation angle varies between 30 and 28 degrees irrespective of the duration of the window Other CVkn criteria than that of the relation (1) Of course, other types of norm than the norm L1 of a vector of a vector can be used, and alternatively a criterion R2 or Nash is used. 16 and 17 are both better as they approach the value of 1, the process seeking this time a maximum and not a minimum as above. at most scarcely marked, which can make his determination more difficult. The criterion according to relation (1), which has a very marked minimum, is therefore preferred.
Figure 18 illustrates a direction of propagation determined manually by an operator. In two images spaced 7 times apart (28 seconds) the operator drew by hand three vectors of cloud movements. The operator has identified in both images two areas corresponding to the same clouds. The origin of a vector thus corresponds to a cloud pixel in the first image and the end of the vector corresponds to the corresponding pixel in the second image. These three measurements of direction of movement of the clouds are illustrated by the three reporters and give values that are substantially identical and coherent with the direction determined by the method according to the invention (30 °), it being understood that the positioning of the vectors by the operator is uncertain to a few pixels because of cloud deformation.
The mean direction of propagation of the clouds Dkoptn is thus known at each instant n, it is possible to implement a determination of their average speed, according to any known treatment of the state of the art. A calculation according to the invention of this speed is now described. C.4. Estimation of the average speed of cloud movement
Once the mean direction of displacement Dkoptn of clouds has been determined at time n, advantageously, but optionally, according to the method described above, the method continues with the estimate 504 of their average speed in said direction and at this moment .
More particularly, this estimation starts with a step 504-2 consisting, for rectified binary image IN £, resulting from the pretreatment of the newly acquired image INn, to determine the distribution DVkopt / n (DV for "velocity distribution") of the pixels of the zone Zkoptn, n (ie the zone of the image IN% associated with the direction Dkoptn) corresponding to the clouds, along the length L of the rectangle Rkopt, n, and thus in the direction Dk.ovt, n · Ce computation is for example carried out on the rectangular binary image lRkopt, n after rotation of the angle 90 ° - akopt by calculating a histogram counting, for each line of the image IRkopt / n after rotation, the number of pixels of the line which are equal to 1, and therefore a histogram in a direction parallel to the direction Dk n, as illustrated in FIG. 19.
A distribution flow FDVkopt§n = [DVkoptil, DVkopti2, -, DVkoptin_ltDVkopt; n ^ is thus obtained by the unit 24, as illustrated in FIG. 20 by adopting the same graphical representation as previously with the abscissa axis. corresponding to the instants n of acquisition. As illustrated in this figure, it can be seen that this stream comprises affine structures that appear intermittently. These structures correspond to cloud passages. They first appear for high ordinates (ie a cloud bank appears at the top of the thumbnail) and then decrease on the y-axis as this cloud bank approaches the solar disk (the center of the solar disk having the coordinate 0), ie decreases in ordinate over successive images in the image stream {INP}. The inventors have found that the slope of these structures is directly related to the average speed V of cloud displacement in the Dkoptin direction (in pixels per second) according to the relation: where β is the angle formed by the direction of the structures affine with the x-axis.
In a subsequent step 504-4, the method therefore determines the value of the angle β of the affine structures as a function of the last q calculated distributions {DVkopt: np, DVkoptin_p + 1,..., DVkopt: n_1, DVkopt ^, ie on a time window [(n - q) x Te; η X Te], where q is a predetermined integer. The duration of this time window is determined, for example during a previous test campaign. For example, a time window between 4 seconds and 44 seconds is chosen, which corresponds to q = 1, 2,..., 11 for an acquisition step Te of 4 seconds, which a priori is sufficient for a determination. with sufficient precision of the average speed, while ensuring that the speed of the clouds does not undergo much modification. Other durations are possible, for example more important. Any suitable algorithm implemented by the unit 24 may be suitable for this purpose (eg a linear regression).
Note however in Figure 20 that the base structures (areas near the coordinate 0) may be more or less wide because of the proximity to the sun. This flare can disturb the determination of the angle by conventional algorithms. In a preferred manner, an approach similar to that used for the identification of the Dkoptn direction is implemented. More particularly, different rotations of the distribution flow FDVk n are tested (corresponding graphically to different rotations of the axes of FIG. 20). Normally, for a certain rotation of these axes, the affine structures become vertical. In this case, the distribution on the abscissa axis of the cloud pixels has peaks of values greater than the peaks of the distributions on the x-axis of the cloud pixels of different rotations (the peaks are the weakest when the axes are rotated a value leading to affine structures parallel to the abscissa axis). For this purpose, the method according to the invention comprises, for example: the production of a matrix DVn = [DVkoptn_q DVkoptq + 1 ... Wfcopt; "_WWcct;"], the DVkoptJl distributions being coded in the form of a column vector of dimension L the application of different rotations of angles βη predetermined and stored in the unit 24, and center for example coordinates (nq, L) to the matrix DVn. The angles βη are angularly spaced between 0 ° and 90 °, with an angular pitch that adjusts the desired accuracy on the desired average speed, for example a step equal to 0.5 °; calculating a Histnu histogram for each matrix DVnu obtained by rotation and calculating the maximum Mnu of this histogram, namely the maximum of the vector [Σ) = ι dvu (l, j) Σ} = ι dvu (2, j ) ... Σ) = ι dvu (nq, j)]; the determination of the angle βηορ (, η as being the angle βη corresponding to the maximum maximum Mnudes Mn u calculated.
The method then calculates at 504-6 the average velocity Vn of the clouds, in pixels per second, based on the relation (2).
The upper portion of FIG. 2f graphically illustrates the rotation of the matrix DVn illustrated in FIG. 20 for an angle βη. FIG. 22 illustrates the evolution of the maximum Mnu as a function of the angle βη. It is found that the maximum of the curve Mnu (/? U) has a marked maximum corresponding to the angle βηορί, η · Figure 23 illustrates the manual determination by an operator of the angle βηορΙ, η which traces the slopes of the principal FIG. 20. The protractor on the right of FIG. 23 shows that the angle obtained by the method according to the invention is consistent with a measurement of it made directly on the image of the distributions DVk n.
In addition, the deviation of the angle /?Uop(.n relative to a reference angle (eg 180 ° if the angles are measured between 0 ° and 360 ° or 0 ° if the angles are measured between -180 ° and 180 °) gives the direction of propagation of the clouds in the direction Dkopt / n C.5.
Once the average speed of the clouds has been determined at the moment n, the process continues, in 506, by the estimate of the sunshine, in particular the estimate of the fall of sunshine due to clouds passing in the cone solar.
Referring to FIG. 24 which illustrates a portion of the FDVk stream n, the ordinate axis represents the distance (in pixels) at the base of the detection rectangles (^ / copt, n), i.e. the distance (in pixels) at the center of the solar disk in the image flow {INn} This distance can be transformed into time since the average moving speed Vn of the clouds is known, for example, by taking again the speed calculated according to of Figure 20 (1.48 pixel / second), a displacement of 100 pixels on the ordinate axis corresponds to a duration of about 100 / 1.48 ~ 68 seconds (about 17 images in the particular case of an acquisition period Te of 4 seconds.) In other words, a cloud detection performed at the ordinate equal to 100 suggests that the cloud will intercept the solar disk 68 seconds later (or 17 images later). in the case of this sequence, the maximum time anticipation is 500 / 1.48 ~ 340 sec waves, ie 5 minutes 40 seconds (500 being the size of the rectified images and the length Z of the detection rectangles).
Still with reference to FIG. 24, it can be seen that the clouds being progressed towards the solar disk correspond to the significantly non-zero values of the DVk n distribution of the last acquired image INn. Knowing the position in pixels of the clouds at this time (ie the ordinate of the non-zero values of the distribution DVkoptlnl ü is thus possible to predict their arrival time in the solar cone and on the solar disk. The invention therefore estimates the passage of clouds in front of the sun as a function of the ordinates of non-zero values of the last distribution DVkopt n and the average speed Vn of cloud displacement.
For example, the estimation of sunshine consists in detecting, at 506-2, the significantly non-zero values of the last distribution DVk n and storing their respective ordinates {Ys}, then, at each of these ordinates Ys, at determine, in 506-4,
Y the instant of arrival Ts = - corresponding clouds in the solar halo. Moments
Vn calculated {Ts} thus correspond to the moments for which it is predicted a fall of Γsoleillementillement. The instants {Ts} are for example transmitted by the unit 24 to a control unit of the photovoltaic park so that it controls accordingly said park. Another use of the arrival times is to use them for the diagnosis of the proper operation of the photovoltaic park. In particular, a park is monitored and alarms are triggered when a drop in production with unknown origins is detected. Knowing the arrival times of clouds, this alarm can be deactivated or classified with a lower alarm level (eg non-critical).
To illustrate the predictive capacity of the process, the chronicles of the number of cloud pixels per fixed ordinate (ie corresponding to a temporal arrival time at the fixed sun position) were extracted from the image of FIG. ), namely 50, 100, 150, 250, 300 and 350 pixels, as illustrated respectively in FIGS. 25A-G. The speed Vn being 1.48 pixels per second, or 5.67 pixels per time step Te, temporally shifting 25A-G correspondingly (by 50 / 5.67 = 8.81 time steps for the Figure 25A, 100 / 5.67 = 17.6 time steps for Figure 25B, 150 / 5.67 = 26.5 time steps for Figure 25C, 200 / 5.67 = 35.3 no time for Figure 25D, 250 / 5.67 = 44.0 time steps for Figure 25E, 300 / 5.67 = 52.9 time steps for Figure 25F, and 350 / 5.67 = 61.7 time steps for Figure 25G), the peaks of the figures thus shifted must match and also correspond to peaks of sunshine drop. Figure 26 corresponds to the superposition of the offset figures. As can be seen, the peaks of the different offset figures are superimposed, further verifying that the speed of the clouds has remained substantially constant. Figure 27 illustrates a signal recorded by a pyrometer disposed near the system for validation of the method according to the invention during the recording of the images having occurred in Figure 20 as well as the superposition of the offset figures. As is known per se, a pyranometer makes it possible to measure the amount of solar energy incident on it. As can be seen, the signal falls of the pyranometer, which correspond to falls in sunshine, coincide with the peaks of the offset figures, thus proving the ability of the process to predict sunshine, and in particular the falls of sunshine .
D. CALIBRATION OF THE ACTUAL POSITION OF THE IMAGING DEVICE
As described above, the optimal positioning of the camera is difficult to obtain, so that it has a real position (0 [, 0j, i4zr) different from the desired theoretical position (θ1, θ2, Αζ). It will now be described the calibration step 402 implemented by the unit 24, whose purpose is to determine the parameters of this correction. This calibration can be performed daily, automatically and / or on command of a user, after movement of the camera for example (following a particular gust). This calibration also makes it possible to detect if the device has moved. In particular, by implementing the calibration and comparing the actual position obtained with a real position obtained during a previous calibration, it is determined that the apparatus has moved if the difference between the two positions is greater than a threshold. If there is a significant difference, an alarm can be triggered to obtain the movement of an operator in order to reposition the camera.
In a first step 402-2, a digital calibration image stream {IN £} is acquired by the system 10 while the sky is clear and the sun is present in the portion of clear sky photographed.
In a next step 402-4, the actual position P "of the center of the solar disk in each INfi image is calculated according to the image itself. In particular, the calculation by the unit 24 of the Pfi position of the center of the solar disc implements a fast processing based on the specificity of the acquired image, namely that of the sky, and thus based on the presence of zones with very strong blue dominance (clear sky) and areas with more gray tones (cloud and sun). For each image of the stream {IN%}, the computation firstly comprises the computation of the brightness matrix ILcn and of the matrix INn'R ~ B = | R (IN £) - B (IN £) \. The latter makes it possible to separate the clear sky from the rest of the image in the first order, the strong differences corresponding in fact to the parts of clear sky, while the slight differences belong to the clouds and the solar disk. A pixel-by-pixel scan of the matrices Monm 0 IN ^ R ~ B and Monm 0 ILcn is then performed and: if the brightness of the pixel is sufficiently large (ie tZ £ (a, b)> su where su is a threshold close to the maximum max {Monm 0 ILcn) of the matrix Monm 0 ILcn, with sü for example between 80% and 90% of max (Monm 0 ILcn) or for example equal to 250, and that the absolute value of rb of the pixel is sufficiently small , meaning that the pixel is grayer than blue (ie ίη £ 'β-β (α, b) <sr_b where sr_b is a threshold close to zero, with sr_b less than or equal to 5), the pixel (a, b) is declared as a candidate for the solar disk. The candidate pixels thus constitute a subset of the image Monm 0 ILcn; This subset is eroded to eliminate the candidate pixel groups of small sizes and to promote the subsequent identification of distinct areas of candidate pixels. A circular erosion filter with a diameter of seven pixels is for example implemented. There is thus obtained a set of candidate pixel areas; the center of gravity and the distance from the center of gravity to a predetermined position a priori (for example the position of the center of the solar disk calculated in the preceding image or the theoretical position Pf of the center of the solar disk calculated by the optical astronomical model in particular ) are then calculated for each candidate pixel area; the zone whose center of gravity is closest to the predetermined position is retained as corresponding to the solar disk; and the largest circle inscribed in the selected area is then calculated, and corresponds to the solar disk in the image. The actual final position of the center of the sun Pf retained for subsequent calculations implemented by the process is then the center of this circle.
FIG. 28 illustrates an example of actual positions Pf determined for the calibration flow {INf}.
The process continues, at 402-6, for each image of the flux {INf}, by calculating the theoretical position Pf (01, 02, Az) of the center of the solar disk calculated by the optical astronomical model, the latter is calculated using the optical astronomical model knowing the timestamp of the images and arbitrarily fixing the values of the three angles (θνθ2, Αζ). By testing a sufficient number of combinations of values of these three angles (θ ±, θ2, Αζ), it is then possible to determine the combination that produces the best comparison between the calculated Pf coordinates of the center of the solar disk by this analytical method and the coordinates Pf of the center of the solar disk by calculated during step 402-4. For example, the value (θ [, 02, Azr] of the angles finally adopted is that minimizing the sum of the absolute differences between the set of coordinates {Pf} and {Pf}. In general, any optimization method, eg implemented by computer, to solve an optimization problem according to the following relation is suitable:
where Ω is a search set, for example a hypercube centered on the position desired by the operator during the installation of the imaging device (eg (0,0,0)) or M3, Vi = 1, .. ., n, Pf is an estimator of the position of the solar disk in the INf image based on an INf image analysis (eg the Pf position described above), and || || is a cost function quantifying a difference between two M3 vectors (eg the norm L1, L2, infinite, etc.).
Figure 29 illustrates the result of this calibration by indicating with an arrow in the image the orientation of the North azimuth (also represented on the compass at the bottom left), as well as a spirit level symbolizing a point of interest. distance to the horizontality of the device. In this figure, the curve in the background corresponds to the positions of the sun Pn and the curve in the foreground corresponds to the positions of the sun P £ with the three angles equal to (0 [, θ ^, Azr). It can be seen that the two angles defining the horizontal deviations of the fish-eye are less than 10 °, ie angles that are indeed very small, given the precision resulting from the experiment (camera mounted on a foot resting on a gravel coating). As we observe, the theoretical and calculated positions previously empirically coincide very well. This way of proceeding can be considered much more precise than attempts to position by hand, for example a photographic camera, in the field and finally simpler and more precise than the conventional method usually based on digital images of regular checkers to stall an empirical distortion model coupled with a perspective projection.
It has been described a preferred embodiment in which the direction and the average speed of the clouds are calculated for each image acquired, or equivalent to each moment of acquisition. As a variant, a stream of digital images is acquired over a predetermined period deemed appropriate for the calculation (eg of a duration of between 4 seconds and 44 seconds with a minimum number of images, for example between 2 and 11) and stored in the processing unit. Once the last image of the flow acquired, the unit calculates the direction and the speed as previously described (division in zones, calculation of distributions, analysis of the variability on the acquired flux, etc ...) according to all the distributions calculated on this entire stream. The method according to the invention is for example implemented regularly (for example every 5 minutes or at each acquired time window) or at the request of a user. This allows for example to release the resources of the processing unit for other types of tasks.
Likewise, additional treatments can be applied. For example, a detection of outliers may be implemented (values induced for example by the passage of a bird, a sheet, etc ...) and the images or the corresponding distributions discarded. The calculation is then performed for example over a larger time window to keep the same number of images or on the same time window but with a smaller number of distributions.
Claims (17)
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
FR1651122A FR3047829B1 (en) | 2016-02-12 | 2016-02-12 | METHOD OF ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE |
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
FR1651122 | 2016-02-12 | ||
FR1651122A FR3047829B1 (en) | 2016-02-12 | 2016-02-12 | METHOD OF ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE |
Publications (2)
Publication Number | Publication Date |
---|---|
FR3047829A1 true FR3047829A1 (en) | 2017-08-18 |
FR3047829B1 FR3047829B1 (en) | 2019-05-10 |
Family
ID=56101571
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
FR1651122A Active FR3047829B1 (en) | 2016-02-12 | 2016-02-12 | METHOD OF ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE |
Country Status (1)
Country | Link |
---|---|
FR (1) | FR3047829B1 (en) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111695529A (en) * | 2020-06-15 | 2020-09-22 | 北京师范大学 | X-ray source detection method based on human skeleton key point detection algorithm |
CN112785678A (en) * | 2019-11-09 | 2021-05-11 | 李本彦 | Sunshine analysis method and system based on three-dimensional simulation |
CN115170619A (en) * | 2022-06-10 | 2022-10-11 | 山东电力建设第三工程有限公司 | Cloud occlusion prediction method based on dense optical flow method |
-
2016
- 2016-02-12 FR FR1651122A patent/FR3047829B1/en active Active
Non-Patent Citations (6)
Title |
---|
JEAN-FRANCOIS LALONDE ET AL: "What Do the Sun and the Sky Tell Us About the Camera?", INTERNATIONAL JOURNAL OF COMPUTER VISION, KLUWER ACADEMIC PUBLISHERS, BO, vol. 88, no. 1, 15 September 2009 (2009-09-15), pages 24 - 51, XP019796007, ISSN: 1573-1405 * |
L. ROUPIOZ ET AL: "Quantifying the impact of cloud cover on ground radiation flux measurements using hemispherical images", INTERNATIONAL JOURNAL OF REMOTE SENSING, vol. 36, no. 19-20, 18 October 2015 (2015-10-18), GB, pages 5087 - 5104, XP055310935, ISSN: 0143-1161, DOI: 10.1080/01431161.2015.1084440 * |
LIANG WANG ET AL: "Fisheye-Lens-Based Visual Sun Compass for Perception of Spatial Orientation", MATHEMATICAL PROBLEMS IN ENGINEERING, vol. 2011, 1 January 2012 (2012-01-01), CH, pages 1 - 15, XP055310716, ISSN: 1024-123X, DOI: 10.1109/34.765658 * |
REDA I ET AL: "Solar position algorithm for solar radiation applications", SOLAR ENERGY, PERGAMON PRESS. OXFORD, GB, vol. 76, no. 5, 1 January 2004 (2004-01-01), pages 577 - 589, XP004492837, ISSN: 0038-092X, DOI: 10.1016/J.SOLENER.2003.12.003 * |
S. QUESADA-RUIZ ET AL: "Cloud-tracking methodology for intra-hour DNI forecasting", SOLAR ENERGY., vol. 102, 1 April 2014 (2014-04-01), GB, pages 267 - 275, XP055307523, ISSN: 0038-092X, DOI: 10.1016/j.solener.2014.01.030 * |
XI-HUI WANG ET AL: "Research of an omnibearing sun locating method with fisheye picture based ontransform domain algorithm", LECTURE NOTES IN CONTROL AND INFORMATION SCIENCES, BERLIN, DE, vol. 345, 16 August 2006 (2006-08-16), pages 1169 - 1174, XP009095097, ISSN: 0170-8643 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112785678A (en) * | 2019-11-09 | 2021-05-11 | 李本彦 | Sunshine analysis method and system based on three-dimensional simulation |
CN112785678B (en) * | 2019-11-09 | 2023-12-12 | 李本彦 | Sunlight analysis method and system based on three-dimensional simulation |
CN111695529A (en) * | 2020-06-15 | 2020-09-22 | 北京师范大学 | X-ray source detection method based on human skeleton key point detection algorithm |
CN111695529B (en) * | 2020-06-15 | 2023-04-25 | 北京师范大学 | X-ray source detection method based on human skeleton key point detection algorithm |
CN115170619A (en) * | 2022-06-10 | 2022-10-11 | 山东电力建设第三工程有限公司 | Cloud occlusion prediction method based on dense optical flow method |
CN115170619B (en) * | 2022-06-10 | 2023-08-15 | 山东电力建设第三工程有限公司 | Cloud shielding prediction method based on dense optical flow method |
Also Published As
Publication number | Publication date |
---|---|
FR3047829B1 (en) | 2019-05-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US10547786B2 (en) | Image processing for turbulence compensation | |
EP3278301B1 (en) | Method of determining a direction of an object on the basis of an image of the object | |
Jiang et al. | Advancing opportunistic sensing in hydrology: A novel approach to measuring rainfall with ordinary surveillance cameras | |
EP3198311B1 (en) | Detection unit and method for identifying and monitoring clouds in an observed area of the sky | |
FR2882160A1 (en) | Video image capturing method for e.g. portable telephone, involves measuring local movements on edges of images and estimating parameters of global movement model in utilizing result of measurement of local movements | |
WO2008099092A2 (en) | Device and method for watching real-time augmented reality | |
FR3047829A1 (en) | METHOD FOR ESTIMATING THE POSITION OF THE SOLAR DISK IN A SKY IMAGE | |
Filhol et al. | Time‐Lapse Photogrammetry of Distributed Snow Depth During Snowmelt | |
EP3200153A1 (en) | Method for detecting targets on the ground and in motion, in a video stream acquired with an airborne camera | |
EP2909671B1 (en) | Method for designing a single-path imager able to estimate the depth of field | |
EP3682199A1 (en) | System for measuring components of solar radiation | |
EP3473000B1 (en) | Image capture method and system using a virtual sensor | |
FR2993689A1 (en) | METHOD FOR DETERMINING PARAMETERS FOR OBSERVING A WATER SURFACE BY IMAGING | |
FR3047830A1 (en) | METHOD FOR DETERMINING THE DIRECTION OF MOVING OBJECTS IN A SCENE | |
FR3047831A1 (en) | METHOD FOR DETERMINING THE SPEED OF MOVING OBJECTS IN A SCENE | |
EP3216213A1 (en) | Method for detecting defective pixels | |
FR3033913A1 (en) | METHOD AND SYSTEM FOR RECOGNIZING OBJECTS BY ANALYZING DIGITAL IMAGE SIGNALS OF A SCENE | |
EP3170205B1 (en) | Device for detecting movement | |
EP3380947B1 (en) | Method and device for forecasting cloudiness by statistical processing of data selected by spatial analysis | |
FR2981149A1 (en) | Aircraft, has attitude measurement device including optical sensor that captures images of stars, where attitude measurement device measures attitude of aircraft at both day and night from images taken by sensor | |
FR3069681A1 (en) | METHOD AND DEVICE FOR IMAGE CAPTURE OF IRIS | |
FR3049066A1 (en) | SYSTEM FOR MONITORING AND DETECTING AN EVENT ON THE GROUND SURFACE BY A SATELLITE CONSTELLATION | |
WO2016203134A1 (en) | Method for estimating and forecasting an energy production indicator for a solar system | |
EP4217694A1 (en) | Infrared thermographic system and associated method | |
CA3113163A1 (en) | Adaptive image acquisition method, system and computer program |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PLFP | Fee payment |
Year of fee payment: 2 |
|
PLSC | Publication of the preliminary search report |
Effective date: 20170818 |
|
AU | Other action affecting the ownership or exploitation of an industrial property right |
Effective date: 20170822 |
|
PLFP | Fee payment |
Year of fee payment: 3 |
|
PLFP | Fee payment |
Year of fee payment: 5 |
|
PLFP | Fee payment |
Year of fee payment: 6 |
|
PLFP | Fee payment |
Year of fee payment: 7 |
|
PLFP | Fee payment |
Year of fee payment: 8 |
|
PLFP | Fee payment |
Year of fee payment: 9 |