FR 3047829 A1 20170818 FR 1651122 A 20160212 1 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 +1h à +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 3047829 2 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 5 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 10 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 15 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. 20 Toutefois, cette approche pose des problèmes tant au niveau du matériel mis en oeuvre 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. 25 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 30 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 35 moteur pour entrainer 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. 3047829 3 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 5 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 10 comme exempte d'objet en déplacement. Ceci nécessite une calibration très fine du traitement et/ou une post-analyse complexe et couteuse 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é. 15 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 20 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 oeuvre du traitement. 25 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 30 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 35 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. 3047829 4 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 5 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 (0f, t51' , Azr) du dispositif d'imagerie, définie par la position (Of , 0) du capteur numérique par rapport à un plan horizontal et par la position d'un axe du capteur par rapport à un azimut (Azr). 10 Selon l'invention, le procédé comprend une étape de détermination de la position prédéterminée (Of , t51' , Azr) du dispositif d'imagerie, comportant : - l'acquisition d'un flux d'images numériques de calibration horodatées {Ig} comprenant chacune le disque solaire ; 15 - pour chaque image In dudit flux {/N79 , le calcul d'une position Prf du centre du disque solaire dans l'image In en analysant l'image In ; - le calcul de la position du dispositif d'imagerie (Of , t51' , Azr) comme solution d'un problème d'optimisation selon la relation : (Of , 02, Azr) = argmin lpf.\ (Pis (01,02,Az) \131 / 4ns. (01, 02, Az)) (01,02,Az)ef2 20 où SZ est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart entre deux vecteurs de 3 Selon un mode de réalisation, les images sont en couleur, et en ce que le calcul d'une 25 position Prf du centre du disque solaire dans l'image In comprend : - le calcul d'une matrice de luminosité I Lci, d'une matrice /Ni,c'R-13 = IRUN7,9 - B(IN7,9 I égale à la valeur absolue de la différence entre des canaux bleu et rouge de l'image In ; et - le calcul de la position Prf dans la matrice I Lcn. 30 En particulier, le calcul d'une position Prf du centre du disque solaire dans l'image In comprend : - le calcul de zones de pixels de la matrice /Ni,c'R-Bcandidats au disque solaire, un pixel (a, b) étant jugé candidat si sa valeur 172c (a, b) dans la matrice /Lci, est 35c,R-B supérieure à un premier seuil prédéterminé et si sa valeur ini, (a, b) dans la matrice /Ni,c'R-13 est inférieure à un second seuil prédéterminé ; et 3047829 5 - 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 5 de déplacements d'objets en mouvement dans un flux d'images numériques {/Nn} d'une scène comprenant du ciel, les images numériques dudit flux étant munies d'un axe de référence Are f, l'estimation comprenant : l'acquisition du flux d'images numériques {/Nn} ; pour chaque image Mn du flux d'images numériques {MO 10 o le calcul d'une image segmentée Me, identifiant les objets en mouvement dans l'image Mn; o le calcul d'un ensemble de zones d'image {Zka,} 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 Are f 15 o pour chaque zone Z km, dudit ensemble de zones {Zka,} , le calcul d'une distribution DDka, des valeurs des pixels de ladite image segmentée Me, compris dans ladite zone Z kat, la distribution DDka, étant calculée le long d'une direction normale à la direction d'analyse Dk de ladite zone Z kat; pour chaque direction d'analyse Dk, le calcul d'un critère de variabilité CVka, de 20 distributions IDDkm_p,DDka,_p+i, ...,DDka,_1,DDk,721 calculées pour des zones Z associées à ladite direction Dk; et la détermination de la direction moyenne Dkoptm de déplacement des objets dans le flux d'images numériques {/Nn} comme étant égale à la direction d'analyse pour laquelle le critère de variabilité CVka, est optimum. 25 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. Il est par ailleurs robuste car il filtre 30 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 oeuvre 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 35 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 3047829 6 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 {In} comprend : 5 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 DDk,', le redressement des images les images numériques {MO . 10 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. 15 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 {MO , et le procédé comprend le calcul de l'ensemble de zones d'image {Zka,} pour chaque image Mn du flux d'images numériques {MO , 20 ledit calcul comprenant : le calcul d'un ensemble de rectangles {Rk,'} 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 {Zka,} comme étant égales respectivement à l'intersection des 25 rectangles {Rka,} avec l'image Mn. 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. 30 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. 35 3047829 7 Avantageusement, le petit côté de chaque rectangle Rkal de l'ensemble de rectangles {Rka,} 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 5 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 10 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é CVka, pour la direction 15 Dk est calculé en fonction des distributions IDDkm_p,DDk,72_p+1, DDk,72_1, DDk,721 associées aux p > 2 dernières images numériques acquises. Selon un mode de réalisation, le critère de variabilité CVka, pour la direction Dk est calculé selon la relation : 20 = Dpk,j11 j=n-p+1 expression dans laquelle j référence la position d'une image numérique dans le flux d'images numériques {/Nn} , n et n-p référence respectivement la nième image et la (np)lème image dans le flux d'images numériques {/Nn} , avec p strictement inférieur à n, k 25 référençant la klème direction Dk, et 11V II est une norme du vecteur V de l'espace l'Zi I étant la dimension des distributions DDk,n. Les inventeurs ont en effet noté que le critère CVka, présente un minimum global en fonction de l'indice k et que ce minimum global correspond à la direction moyenne de 30 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. 3047829 8 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 DDk,', ledit masque est appliqué aux images du flux d'images numériques {/Nn} et/ou aux images segmentées {/Nin . En d'autres termes, en excluant les objets de la scène qui ne 5 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 10 numériques {/Nn} 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 Dkoptl. Notamment, les objets en mouvement sont des nuages, et le procédé comprend la 15 prévision d'ensoleillement en fonction de nuages dans le flux d'images numériques {/Nn} 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. 20 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- 25 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 /Ne., le calcul d'une distribution DVkopt,n des valeurs 30 des pixels de l'image segmentée Me, appartenant à la zone Zkoptm associée à la direction moyenne déterminée Dkoptal, la distribution DVkopt,n étant calculée le long de la direction moyenne Dkoptal déterminée ; le calcul de l'inclinaison d'au moins une structure affine comprise dans une matrice DVn = [Dilkopt,n_qpilkopt,n-q+1 --- DVkopt,n_1DVkopt,n] issue de la concaténation des 35 distributions DVkopt,n des images entre la (n-q)lème image du flux jusqu'à la nième du flux d'images numériques {/Nn} , 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. 3047829 9 En d'autres termes, les inventeurs ont constaté qu'il apparait naturellement des structures affines dans la matrice DV, = [DVkopt,72_,IDVkopt,'_q+1 ...DVkopt,72_1DVkopt,72]. 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 5 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 DV, comprend : 10 l'application de différentes rotations d'angles le' prédéterminés à la matrice DVn; le calcul d'un histogramme Histim, pour chaque matrice DV,,' obtenue par rotation et le calcul du maximum M72,' de cet histogramme ; la détermination d'inclinaison comme étant égal à l'angle &opta, correspondant au maximum des maxima M72,' calculés. 15 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 20 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 25 (Of , t51' , Azr) du dispositif d'imagerie, définie par la position (Of , 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 (Of , 02, Azr) du dispositif d'imagerie : 30 pour chaque image In d'un flux d'images numériques de calibration horodatées {/N79 comprenant chacune le disque solaire {/N79 , en calcul d'une position Prf du centre du disque solaire dans l'image In en analysant l'image In ; en calculant la position du dispositif d'imagerie (Of , t51' , Azr) comme solution d'un problème d'optimisation selon la relation : 35 (Of , 02, Azr) = argmin pf.\ (Pis (01,02, Az) (01,02,Az)612 \1311`: / Pn (01,02, AZ)) 3047829 10 où D. est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart entre deux vecteurs de 3 En particulier, le système est configuré pour mettre en oeuvre un procédé du type précité. 5 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 10 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 (Of , 0 , Azr) du dispositif d'imagerie, définie par la position (Of , OD du capteur numérique par rapport à un plan 15 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 (Of , 0 , Azr) du dispositif d'imagerie, comportant : - l'acquisition d'un flux d'images numériques de calibration horodatées {Ig} comprenant chacune le disque solaire ; 20 - pour chaque image I n dudit flux {/N79 , le calcul d'une position Prf du centre du disque solaire dans l'image In en analysant l'image In ; - le calcul de la position du dispositif d'imagerie (Of , Oz, Azr) comme solution d'un problème d'optimisation selon la relation : (Of, 02, Azr) = argmin pf. (Pis (01,02,Az) (01,02 ,Az)e 12 \131 / Pn (01, 02, AZ)) 25 où D. est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart entre deux vecteurs de II: 3 En particulier, le support comprend des instructions pour la mise en oeuvre d'un procédé 30 du type précité. 3047829 11 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 : 5 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 oeuvre, par exemple mis en oeuvre par le système de la figure 1 ; 10 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 15 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 ; 20 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 25 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 30 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 ; 35 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 ; 3047829 12 - 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 ; 5 - 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 10 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. 15 - 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. 20 DESCRIPTION DETAILLEE DE L'INVENTION A. NOTATIONS positifs, avec notamment pour des matrices M et N de I':AxB M = (Ma,b)aE [1,A],bE [1,B] MI = (Inta,b DaE[1,A],bE[1,B] 30 La notation M N se rapporte à l'application d'un masque d'exclusion M à 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 à 35 l'objet M N seul (la « matrice M N »), 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 25 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 AxB où A et B sont des entiers I 3047829 13 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 5 masque M de valeur 1. B. MA 1 ERIEL En se référant à la figure 1, un dispositif 10 d'imagerie et de traitement d'image selon 10 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 15 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 20 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 25 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. 30 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 35 oeuvre par ordinateur selon l'invention. 3047829 14 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 oeuvre l'invention. L'unité 24 est 5 notamment pourvue de l'ensemble des mémoires (RAM, ROM, cache, mémoire de masse,...) pour la mémorisation des instructions mettant en oeuvre le procédé selon l'invention, des paramètres et réglages utiles à cette mise en oeuvre et pour la mémorisation des résultats des calculs intermédiaires et finaux. 10 Optionnellement, l'appareil photographique 12 est avantageusement logé dans un boitier 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 boitier de protection 26 et communique, via par exemple un réseau sans fil (téléphonique, satellitaire, etc...), le 15 résultat du procédé qu'elle met en oeuvre à un site distant où la production du parc photovoltaïque est gérée. En fonctionnement, un flux d'images numériques horodatées {MO est donc produit par l'appareil 12, une image numérique IN 7, étant acquise par exemple toutes les Te secondes, 20 l'entier n représentant le nième instant d'acquisition n x Te depuis un instant initial n = 1, et donc la nième image dans le flux d'images {MO . 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 IN?, comprend A lignes de pixels par B colonnes de pixels et comprend 25 trois canaux d'un espace colorimétrique, par exemple les canaux rouge (r), vert (y) et bleu (b) codés respectivement par trois matrices R, V, B deAxB 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). 30 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, 35 l'invention s'applique à un nombre de bits quelconque. 3047829 15 C. TRAI l'EMENT 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 5 photovoltaïque, et plus particulièrement la chute d'ensoleillement sur le parc photovoltaïque entrainant une chute de la production, mis en oeuvre 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. 10 C.1. Calibration Le procédé débute par une phase 40 de positionnement et de calibration du dispositif d'imagerie 10. 15 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 20 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 « Geometric modelling and calibration of fisheye lens camera system », Schwalbe E., In Proc. 2nd 25 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. 30 En connaissant la position (01,612, Az) 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 35 « modèle astronomique optique ». 3047829 16 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, 02,Az), le couple d'angles (01, 02) 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 5 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 à (01, 02,Az) = (0,0,0). Les coordonnées géographiques, e.g. la longitude et la latitude, du dispositif sont par ailleurs 10 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 15 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é. 20 Une fois le dispositif 10 positionné, une calibration 404 est mise en oeuvre. Notamment, le positionnement parfait de l'appareil est difficile à obtenir, de sorte qu'il présente une position réelle (Of , 0 , Azr) différente de la position théorique (01, 02, Az). 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 oeuvre. L'étape de calibration 402 25 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 Of , t51' , 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. 30 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 oeuvre toute méthode connue permettant de déterminer ce type de masque, y 35 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. 3047829 17 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 5 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 IN,,., en journée, par exemple 10 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é I de l'image IN,,.,, d'une manière connue en soi, par exemple selon la formule IL,,, = -31 x + V + B), obtenant ainsi une image en niveau de gris telle qu'illustrée à la figure 3 ; 15 c) la détermination d'un premier masque M1 E 11AxB correspondant à l'image optique de la scène sur le capteur, à savoir mi (a, b) = 1 si le pixel (a, b) de IN,,., 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 selection method form gray-level 20 histograms », de Otsu N, Automatica, 11(285-296), 23-27, 1975. Ce type d'algorithme consiste à diviser l'histogramme de la matrice M1 11,,, 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: (S1, S2, S3) = 25 argmin (6(0, s1) + o-(s1 + 1, s2) + a(s2 + 1, s3) + a(s3 + 1,255)) (sl,s2,S3)E [0,255] 3 où o-(n1, n2) est l'écart type des valeurs de l'histogramme d'abscisses comprises entre n1 et n2. En adoptant 5 niveaux de gris NGI=0 1 pas de temps, à savoir un critère de dissemblance CVk,n pour m = 1 selon la relation : CVk,n = - ddk,i(i))1 (1) j=n-p+1i=1 Les inventeurs ont découvert qu'un tel critère comprend deux optimums (deux maxima si 10 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 15 ak,n 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. 20 minimise) le critère CVk,n entre 0° et 180°, ainsi que l'indice k', correspondant à l'angle ak + 180°, qui maximise (resp. minimise) le critère CVk,n entre 180° et 360°. Par exemple, pour un pas d'angle 61 de 1°, le procédé détermine l'indice k=kopt compris entre 0 et 180 qui maximise (resp. minimise) le critère CVk,n = CVk,n + CVk+180,n. La direction Dkopt,n correspondant à l'indice kw,t déterminé par l'unité 24 pour l'instant n, est ainsi la 25 direction moyenne de déplacement des nuages à l'instant n. La figure 15 illustre des critères CVk,n se fondant sur la relation (1) pour différentes écartements temporels entre les distributions. Les ronds indiquent la valeur de l'angle akopt pour laquelle le critère CVk,n est minimal, et donc la direction Dkopt,n de 30 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. 35 3047829 25 D'autres critères CVka, que celui de la relation (1) peuvent bien entendu être utilisés. Notamment d'autres types de norme que la norme Ll 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 5 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é. 10 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 15 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. 20 La direction moyenne de propagation des nuages Dkopt,n étant ainsi connue à chaque instant n, il est possible de mettre en oeuvre 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. 25 C.4. Estimation de la vitesse moyenne de déplacement des nuages Une fois déterminée la direction moyenne de déplacement Dkopt,n des nuages à l'instant n, avantageusement, mais optionnellement, selon le procédé décrit ci-avant, le procédé se 30 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 INd, issue du prétraitement de l'image nouvellement acquise INn, à 35 déterminer la distribution DVkopta, (DV pour « distribution de vitesse ») des pixels de la zone Zkopt,n,n (i.e. la zone de l'image INd associée à la direction Dkopt,n) correspondant aux nuages, le long de la longueur L du rectangle Rkopt,n, et donc selon la direction Dkopt,n. Ce calcul est par exemple réalisé sur l'image binaire rectangulaire /Rkopt,n après 3047829 26 rotation de l'angle 90° - akop, en calculant un histogramme comptabilisant, pour chaque ligne de l'image /Rkopt,' 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 Dkoptal, tel qu'illustré à la figure 19. 5 Un flux de distributions FDVkopta, = {DVkopt,i, DVkopt,2, ..., DVkopt,n_i, DVkopt,,,} 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 10 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 15 dans le flux d'images {Min. 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 7, (en pixels par seconde) selon la relation : Dkopt, 1 (2) V = tan /3 20 où /3 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 le des structures affines en fonction des q dernières distributions calculées Vkopt,n_p, DVkopt,n-p+1, --- DVkopt,n-1, DVkoptall, i.e. sur une fenêtre temporelle 25 [(n - q) x Te; n 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, 30 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 oeuvre par l'unité 24 peut convenir à cet effet (par exemple une régression linéaire). 35 3047829 27 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 5 direction Dkoptal est mise en oeuvre. Plus particulièrement, différentes rotations du flux de distribution FDVkopt,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 10 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 : 15 la production d'une matrice DVn = [DVkopt,n_q DVkopt,q+1 ...DVkopt,n-ipVkopt,n1, les distribution DVkopt,n étant codées sous la forme de vecteur colonne de dimension L l'application de différentes rotations d'angles leu 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 leu sont espacés angulairement entre 0° et 90°, avec un pas angulaire qui règle 20 la précision souhaitée sur la vitesse moyenne recherchée, par exemple un pas égal à 0,5° ; le calcul d'un histogramme Histnx pour chaque matrice DVnx obtenue par rotation et le calcul du maximum Mn,' de cet histogramme, à savoir le maximum du vecteur [E=1 dvu(1,i) dvu(2,i) dvu(n q,j)]; 25 la détermination de l'angle 8 uopt,n comme étant l'angle leu correspondant au maximum maxis Mn,udes Mnx 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). 30 La portion supérieure de la figure 21 illustre graphiquement la rotation de la matrice DVn illustrée à la figure 20 pour un angle Pu. La figure 22 illustre l'évolution du maximum Mnx en fonction de l'angle Pu. On constate que le maximum de la courbe Mn,u(flu) présente un maximum marqué correspondant à l'angle fluopt,n. La figure 23 illustre la 35 détermination manuelle par un opérateur de l'angle fluopt,n qui trace les pentes des principales structures de la figure 20. Le rapporteur à droite de la figure 23 montre que l'angle flopt,nobtenu par le procédé selon l'invention est cohérent avec une mesure de celui-ci faite directement sur l'image des distributions DVkopt,n. 3047829 28 En outre, l'écart de l'angle /3 01 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 Dkoptal. 5 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. 10 En se référant à la figure 24 qui illustre une partie du flux FDVkopt,', l'axe des ordonnées représente la distance (en pixels) à la base des rectangles de détection {Rkopt,'} , c'est-à- dire la distance (en pixels) au centre du disque solaire dans le flux d'images {/Nn} . Cette distance peut être transformée en temps puisque la vitesse moyenne de déplacement 15 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 20 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 L des rectangles de détection). 25 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 DVkopt,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 DVkopt,n), il est ainsi possible de prédire leur moment d'arrivée dans le cône solaire et sur 30 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 35 significativement non nulles de la dernière distribution DVkopt,n et à mémoriser leurs ordonnées respectives {Ys} , puis, à chacune de ces ordonnées Ys, à déterminer, en 506-4, Ys l'instant d'arrivée Ts = V- des nuages correspondant dans le halo solaire. Les instants n calculés {Ts} correspondent ainsi aux instants pour lesquels il est prédit une chute de 3047829 29 l'ensoleillement. Les instants {n} 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 5 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 10 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 15 temps pour la Figure 25A, 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 à 20 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 25 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 30 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 35 difficile à obtenir, de sorte qu'il présente une position réelle (Of , t51' , Azr) différente de la position théorique souhaitée (01,612, Az). Il va à présent être décrit l'étape de calibration 402 mis en oeuvre 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 3047829 30 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 oeuvre la calibration et en comparant la position réelle obtenue avec une position réelle obtenue au cours d'une 5 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. 10 Dans une première étape 402-2, un flux d'images numériques de calibration {/A/7" 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 Prf du centre du disque solaire dans 15 chaque image INn est calculée en fonction de l'image elle-même. Notamment, le calcul par l'unité 24 de la position Prf du centre du disque solaire met en oeuvre 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 {/A/7" , le calcul 20 comprend tout d'abord le calcul de la matrice de luminosité aci, et de la matrice nvnc,R-B = IR(INnc) - B (INnc) I. 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 nvnc,R-B et monm balayage pixel par pixel des matrices M'in /Lcn est ensuite 25 effectué et : si la luminosité du pixel est suffisamment importante (i.e. Un'. (a, b) > su où su est un seuil proche du maximum max(monm acn) de la matrice M'in acn, avec su compris par exemple entre 80% et 90% de max( monm I Lcn) ou par exemple égal à 250, et que la valeur absolue de r-b du pixel est suffisamment petite, signifiant que 30 le pixel est plus gris que bleu (i.e. illnc'R-B (a, b) REVENDICATIONS1. 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 (Of , Azr) du dispositif d'imagerie, définie par la position (01, OD du capteur numérique par rapport à un plan horizontal et par la position d'un axe du capteur par rapport à un azimut (Azr), caractérisé en ce que le procédé comprend une étape de détermination de la position prédéterminée (01, Azr) du dispositif d'imagerie, comportant : l'acquisition d'un flux d'images numériques de calibration horodatées UN79 comprenant chacune le disque solaire ; pour chaque image 11^17.c, dudit flux f/N79 , le calcul d'une position Prf du centre du disque solaire dans l'image INrci en analysant l'image INn ; le calcul de la position du dispositif d'imagerie (Of , Azr) comme solution d'un problème d'optimisation selon la relation : (Of , Azr) = argmin (Pi ) (Pis (01, 02, Az) (01,192,Az)ell (91, 02, AZ)) où S1 est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart entre deux vecteurs de I83. 2. Procédé selon la revendication 1, caractérisé en ce que les images sont en couleur, et en ce que le calcul d'une position 13,F du centre du disque solaire dans l'image INfi comprend : le calcul d'une matrice de luminosité 14; d'une matrice /N7ic'R-B = IR UND - B(//\179 I égale à la valeur absolue de la différence entre des canaux bleu et rouge de l'image le, ; et le calcul de la position /Sri dans la matrice I Lcn. 3047829 34 3. Procédé selon la revendication 2, caractérisé en ce que le calcul d'une position Pif du centre du disque solaire dans l'image INn comprend : - le calcul de zones de pixels de la matrice /Nnc'R-Bcandidats au disque solaire, un pixel (a, b) étant jugé candidat si sa valeur ln (a, b) dans la matrice I Ln est 5 supérieure à un premier seuil prédéterminé et si sa valeur inn'R-B (a, b) dans la matrice /Nnc'R-B 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. 10 4. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il 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!, ladite estimation comprenant : 15 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 Min identifiant les objets en mouvement dans l'image INn; o le calcul d'un ensemble de zones d'image tZk,n1 respectivement 20 associées à des directions d'analyse {Dk} formant des angles fak} différents prédéterminés avec l'axe de référence Are f : o pour chaque zone Zki, dudit ensemble de zones {Ziui} , le calcul d'une distribution DDk,n des valeurs des pixels de ladite image segmentée In compris dans ladite zone Zkx, la distribution DDk,i, étant calculée le long 25 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é CVk,n de distributions {DDk,n_p, Dpk,n-p+1, Dpk,n-1, Dpic,n} calculées pour des zones Zk,n associées à ladite direction Dk; et la détermination de la direction moyenne Dkopt,n de déplacement des objets 30 dans le flux d'images numériques {INn} comme étant égale à la direction d'analyse pour laquelle le critère de variabilité CVk,' est optimum. 5. Procédé selon la revendication 4, caractérisé en ce que l'acquisition du flux d'image numériques {INn} 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° ; 3047829 35 - et préalablement au calcul des distributions DDk,', le redressement des images les images numériques {INn} . 6. Procédé selon la revendication 4 ou 5, caractérisé en ce que la scène est acquise à 5 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 en ce que le procédé comprend le calcul de l'ensemble de zones d'image {Zka,} pour chaque image INn du flux d'images numériques {INn} , ledit calcul comprenant : 10 - le calcul d'un ensemble de rectangles {Rk,72} 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 {Zka,} comme étant égales respectivement à l'intersection des rectangles {Rka,} avec l'image INn. 15 7. Procédé selon la revendication 6, caractérisé en ce que la scène comprend du ciel, et en ce que le point prédéterminé du plan focal est la position du centre du disque solaire sur le plan focal du système optique. 20 8. Procédé selon la revendication 7, caractérisé en ce le petit côté de chaque rectangle Rkal de l'ensemble de rectangles {Rk,72} est au moins égal au diamètre du halo D7211 du disque solaire. 9. Procédé selon l'une quelconque des revendications 4 à 8, caractérisé en ce que le 25 critère de de variabilité CVka, pour la direction Dk est calculé en fonction des distributions IDDka'p, DDka'p+1, , DDk,72_1, DDk,721 associées aux p > 2 dernières images numériques acquises. 10. Procédé selon l'une quelconque des revendications 4 à 9, caractérisé en ce que le 30 critère de variabilité CVka, pour la direction Dk est calculé selon la relation : C1,1k,n - IIDDki_i DDkill j=n-p+1 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 nième image et la (n-p)lème image dans le flux d'images numériques {INn} , avec p strictement 35 inférieur à n, k référençant la klème direction Dk, et 11VII est une norme du vecteur V de l'espace l'Zi , I étant la dimension des distributions DDk,n. 3047829 36 11. Procédé selon l'une quelconque des revendications 4 à 10, caractérisé en ce qu'il comporte le calcul d'un masque d'exclusion des éléments immobiles de la scène, et en ce que préalablement au calcul des distributions DDk,', ledit masque est appliqué aux images du flux d'images numériques {/Nn} et/ou aux images 5 segmentées {Min . 12. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il comprend en outre : le calcul d'une vitesse moyenne Vn des objets en mouvement dans le flux 10 d'images numériques {/Nn} 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 Dkoptm. 15 13. Procédé selon la revendication 12, caractérisé en ce qu'il comprend la prévision de l'ensoleillement au-dessus d'un parc photovoltaïque, en particulier la chute d'ensoleillement. 14. Procédé selon la revendication 12 ou 13, caractérisé en ce que le calcul de ladite 20 vitesse moyenne Vn comprend : pour chaque image segmentée INe., le calcul d'une distribution DVkopt,n des valeurs des pixels de l'image segmentée INe, appartenant à la zone Zkoptal associée à la direction moyenne déterminée Dkopt,n, la distribution DVkopt,n étant calculée le long de la direction moyenne Dkoptal déterminée ; 25 le calcul de l'inclinaison d'au moins une structure affine comprise dans une matrice DVn = [DVkopt,n_qDVkopt,n_q+1 DVkopt,n_1DVkopt,n1 issue de la concaténation des distributions DVkopt,n des images entre la (n-q)lème image du flux jusqu'à la nième du flux d'images numériques {/Nn} , où q est un entier supérieur ou égal à 1; 30 et le calcul de la vitesse moyenne Vn en fonction de l'inclinaison calculée. 15. Procédé selon la revendication 14, caractérisé en ce que le calcul de l'inclinaison de la au moins une structure affine comprise dans une matrice DV, comprend : l'application de différentes rotations d'angles leu prédéterminés à la matrice 35 DV, ; le calcul d'un histogramme Histnx pour chaque matrice DV,,' obtenue par rotation et le calcul du maximum Mn,' de cet histogramme ; 3047829 37 - la détermination d'inclinaison comme étant égal à l'angle &opta, correspondant au maximum des maxima M72,' calculés. 16. Système d'estimation de la position du disque solaire dans une image numérique : 5 - 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 (Of , Oz, Az') du dispositif d'imagerie, définie par la position 10 (Of , OD 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 (Of, 0z, Az') du dispositif d'imagerie : - pour chaque image In d'un flux d'images numériques de calibration 15 horodatées {/N79 comprenant chacune le disque solaire {/N79 , en calcul d'une position Prf du centre du disque solaire dans l'image In en analysant l'image In ; - en calculant la position du dispositif d'imagerie (Of, 0z, Azr) comme solution d'un problème d'optimisation selon la relation : 20 (Of , 02, Azr) = argmin pf. (Pis (01,02, Az) (01,02,Az)612 \131 / 4ns. (01, 02, AZ)) où D. est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart entre deux vecteurs de 3 25 17. Système selon la revendication 16 configure pour mettre en oeuvre un procédé selon l'une quelconque des revendications 2 à 15.