FR 3047830 A1 20170818 FR 1651120 A 20160212 1 PROCEDE DE DETERMINATION DE LA DIRECTION DE DEPLACEMENT D'OBJETS DANS UNE SCENE DOMAINE DE L'INVENTION L'invention a trait au domaine du traitement d'un flux d'images numériques d'une scène en vue de déterminer la direction moyenne de déplacement d'objets dans la scène. L'invention trouve particulièrement application dans la détermination de la direction des nuages dans un ciel. La direction ainsi déterminée 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 3047830 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. 3047830 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 EXPOSE DE L'INVENTION Le but de la présente invention est de proposer un procédé d'estimation de la direction moyenne de déplacement d'objets en mouvement dans une scène qui soit robuste aux phénomènes locaux, peu consommatrice de ressources informatiques, et qui peut 30 s'employer le cas échéant également avec une image produite avec une optique à très grand angle de champ, y compris une optique de type « fisheye ». A cet effet, l'invention a pour objet un procédé d'estimation de la direction moyenne de déplacements d'objets en mouvement dans un flux d'images numériques {MO d'une 35 scène comprenant du ciel, les images numériques dudit flux étant munies d'un axe de référence Are f, le procédé comprenant : - l'acquisition du flux d'images numériques {Mn} ; - pour chaque image Mn du flux d'images numériques {MO 3047830 4 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 5 prédéterminés avec l'axe de référence A' f : 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 Zkal; 10 pour chaque direction d'analyse Dk, le calcul d'un critère de variabilité CVka, de 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 15 laquelle le critère de variabilité CVka, 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 20 adapté à la course des nuages dans le ciel. Il 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 oeuvre sont basiques (calcul de distributions et calculs des variations de distribution), et ne nécessitent donc pas de 25 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 30 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 : l'acquisition de la scène par une dispositif de prise de vue ayant un angle de champ 35 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 {/Nn} . 3047830 5 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. 5 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 {/Nn} , 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 {/Nn} , 10 ledit calcul comprenant : 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 15 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. 20 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. 25 Avantageusement, le petit côté de chaque rectangle Rical de l'ensemble de rectangles {Ric,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. 30 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 35 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. 3047830 6 Plus particulièrement, le modèle astronomique est réglé pour une position prédéterminée (Of , 0z, Azr) du dispositif d'imagerie, définie par la position (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 en ce que le procédé comprend une étape de détermination 5 de la position prédéterminée (Of , 0z, 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 ; - 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 ; 10 - le calcul de la position du dispositif d'imagerie (Of , , Azr) comme solution d'un problème d'optimisation selon la relation : (Of, 02 , Azr) = argmin pf.\ (Pis (01,02,Az) (01,02,Az)612 \131 / 4ns. (01, 02, Az)) où SZ est un ensemble de recherche, et II II est une fonction de coût quantifiant un 15 écart entre deux vecteurs de 3 En d'autres termes, malgré tout le soin que l'on peut mettre dans le positionnement du dispositif d'imagerie, il peut exister une erreur entre la position réelle du dispositif et la position souhaitée pour celui-ci. L'invention permet de corriger simplement, 20 automatiquement et avec précision cette erreur de positionnement. Notamment, les images sont en couleur, et en ce que le calcul d'une position Prf du centre du disque solaire dans l'image In comprend : - le calcul d'une matrice de luminosité I Lci, d'une matrice /Ni2c'R-13 = IRUN7.9 - 25 B(I1\17,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. En particulier, le calcul d'une position Prf du centre du disque solaire dans l'image In 30 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 172c (a, b) dans la matrice /Lci, est supérieure à un premier seuil prédéterminé et si sa valeur ini,c'R-B (a, b) dans la matrice /Ni2c'R-13 est inférieure à un second seuil prédéterminé ; et 35 - la sélection d'une zone candidate comme appartenant au disque solaire en fonction de sa proximité à une position prédéterminée. 3047830 7 Selon un mode de réalisation particulier, le critère de de variabilité CVka, pour la direction Dk est calculé en fonction des distributions {DDk,n_p,DDk,n_p+i, DDk,n} associées aux p > 2 dernières images numériques acquises. 5 Selon un mode de réalisation, le critère de variabilité CVka, pour la direction Dk est calculé selon la relation : CVkal = Dpk,j11 j=n-p+1 expression dans laquelle j référence la position d'une image numérique dans le flux 10 d'images numériques {/Nn} , 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 {/Nn} , avec p strictement inférieur à n, k 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. 15 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 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. 20 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,n, 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 25 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 30 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 Dkopt,n. 3047830 8 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 {/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 5 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 10 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. 15 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 des pixels de l'image segmentée Me, appartenant à la zone Zkopt,n associée à la direction moyenne déterminée Dkopt,', la distribution DVkopt,n étant calculée le long de la direction moyenne Dkopt,n déterminée ; 20 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,n] issue de la concaténation des distributions DVkopt,n des images entre la (n-q)ieme 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. 25 En d'autres termes, les inventeurs ont constaté qu'il apparait naturellement des structures affines dans la matrice DVn = [DVkopt,n_qDVkopt,n_q+1 DVkopt,n_iDVkopt,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 30 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 : 35 l'application de différentes rotations d'angles le' prédéterminés à la matrice DVn ; le calcul d'un histogramme Histn,' pour chaque matrice DVn,n obtenue par rotation et le calcul du maximum Mn,n, de cet histogramme ; 3047830 9 - la détermination d'inclinaison comme étant égal à l'angle &opta, correspondant au maximum des maxima M72,' calculés. Ce calcul permet de déterminer de manière précise l'inclinaison à l'aide de calcul 5 algébriques simples, et donc peu consommateurs de ressources informatiques. L'invention a également pour objet un système d'estimation de la direction moyenne de déplacements d'objets en mouvement dans un flux d'images numériques {MO d'une scène comprenant du ciel, les images numériques dudit flux étant munies d'un axe de 10 référence Are f, le système comprenant : une caméra pour l'acquisition du flux d'images numériques {MO de ladite scène; une unité de traitement d'information configurée pour : o pour chaque image Mn du flux d'images numériques {MO ^ le calcul d'une image segmentée Me, identifiant les objets en mouvement 15 dans l'image Mn; ^ 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 A' f ^ pour chaque zone Z km, dudit ensemble de zones {Zka,} , le calcul d'une 20 distribution DDk,i, 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 km; o pour chaque direction d'analyse Dk, le calcul d'un critère de variabilité CVka, de distributions {DDka,_p, DDka,_p+i, , DDk,72_1, DDka,} calculées pour des 25 zones Z km associées à ladite direction Dk; et o la détermination de la direction moyenne Dkoptal 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. 30 En particulier, le système est configure pour mettre en oeuvre 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 direction moyenne de déplacements d'objets 35 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, le procédé comprenant : 3047830 10 pour chaque image Mn du flux d'images numériques {MO : 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 à 5 des directions d'analyse {Dk} formant des angles {ak} différents prédéterminés avec l'axe de référence A' f : 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 10 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 distributions IDDkm_p,DDka,_p+i, DDk,72_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 15 flux d'images numériques {/Nn} comme étant égale à la direction d'analyse pour laquelle le critère de variabilité CVka, est optimum. En particulier, le support comprend des instructions pour la mise en oeuvre d'un procédé du type précité. 20 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 : 25 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 ; 30 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 35 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 ; 3047830 11 - 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 ; 5 - 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 ; 10 - 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 15 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 20 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 ; 25 - 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 30 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. 35 - 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 ; 3047830 12 - 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 5 A. NOTATIONS 10 positifs, avec notamment pour des matrices M et N deAxB M = (Ma,b)aE [1,A],bE [1,B] IM I = (Inta,b DaE[1,A],bE[1,B] La notation MON se rapporte à l'application d'un masque d'exclusion M à la matrice N. Plus particulièrement, l'application d'un traitement sur l'objet MON signifie que 15 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 M ON seul (la « matrice MON»), 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 20 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 » MON 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. 25 B. MA I ERIEL 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 30 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 35 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 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 3047830 13 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 5 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 10 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 15 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 oeuvre par ordinateur selon l'invention. 20 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 notamment pourvue de l'ensemble des mémoires (RAM, ROM, cache, mémoire de 25 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. Optionnellement, l'appareil photographique 12 est avantageusement logé dans un boitier 30 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 résultat du procédé qu'elle met en oeuvre à un site distant où la production du parc 35 photovoltaïque est gérée. 3047830 14 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 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 {INn} . En variante, l'acquisition est réalisée 5 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 (y) et bleu (b) codés respectivement par trois matrices R, V, B de 11AxB Les valeurs numériques 10 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 15 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. 20 C. TRAIIEMENT 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 25 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. C.1. Calibration 30 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 35 instructions et paramètres utilisés par le procédé, notamment : 3047830 15 - 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 5 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 Panaramic Photogrammetry Worshop, Int. Archives of Photogrammetry and Remote 10 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 (01, 02, Az) du capteur 14 dans un repère terrestre, et donc la 15 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 ». 20 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 25 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 30 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 35 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é. 3047830 16 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 ,t51',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 5 l'écart entre les positions réelle et théorique est mise en oeuvre. 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 (Of, 0z, 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. 10 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 15 de mettre en oeuvre 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. 20 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é 25 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 en journée, par exemple par temps ensoleillé ou légèrement nuageux avec le soleil haut dans le ciel, e.g. à son 30 zénith ; b) le calcul d'une matrice de luminosité de l'image IN,' d'une manière connue en soi, par exemple selon la formule = -31 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 M1 E 11AxB correspondant à l'image optique 35 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 ; 3047830 17 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 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 5 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) = argmin (6(0, s1) + o-(s1 + 1, s2) + a(s2 + 1, s3) + a(s3 + 1,255)) (sl,s2,53) E [0,255] 3 10 où a(nl,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 20 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. 25 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 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. 30 Le procédé détermine ainsi l'indice k, correspondant à l'angle ak, qui maximise (resp. 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 35 0 et 180 qui maximise (resp. minimise) le critère CVk,n = CVk,n + CVk+180,n. La direction 3047830 24 Dkopt,n correspondant à l'indice If0pt 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 CVka, se fondant sur la relation (1) pour différentes 5 écartements temporels entre les distributions. Les ronds indiquent la valeur de l'angle akopt pour laquelle le critère CVkal est minimal, et donc la direction Dkoptal 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 10 28 degrés quel que soit la durée de la fenêtre. 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 15 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é. 20 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 25 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 30 incertain à quelques pixels près en raison de la déformation des nuages. La direction moyenne de propagation des nuages Dkopt,p é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 35 vitesse est à présent décrit. 3047830 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 5 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 INe., issue du prétraitement de l'image nouvellement acquise INn, à 10 déterminer la distribution DVkopta, (DV pour « distribution de vitesse ») des pixels de la zone Z kopt n (i.e. la zone de l'image INe, associée à la direction D kopt,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 /Rkopta, après rotation de l'angle 90° - akopt en calculant un histogramme comptabilisant, pour chaque 15 ligne de l'image /Rkopt,i, 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 Dkopt,n, tel qu'illustré à la figure 19. Un flux de distributions FDVkopta, = {DVkopt,i,DVkoptDVkopt,n} est ainsi 20 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 25 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 {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 30 Dk"t,n (en pixels par seconde) selon la relation : 1 V = tan /3 (2) le est l'angle que forme la direction des structures affines avec l'axe des abscisses. 35 3047830 26 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+i, --- ,DVkopt,n-1, DVkopt,n}, i.e. sur une fenêtre temporelle [(n - q) x Te; n x Te], où q est entier prédéterminé. La durée de cette fenêtre temporelle 5 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. 10 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). On note cependant sur la figure 20 que la base des structures (zones proches de la 15 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 Dkopt,n est mise en oeuvre. Plus particulièrement, différentes rotations du flux de distribution FDVkopt,n sont testées (correspondant graphiquement à différentes rotations 20 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 25 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 = [DVkopt,n_q DVkopt,q+1 ...DVkopt,n-POVicept,n1, les distribution DVkopt,n étant codées sous la forme de vecteur colonne de dimension L 30 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 la précision souhaitée sur la vitesse moyenne recherchée, par exemple un pas égal à 05°; 35 le calcul d'un histogramme Histnx pour chaque matrice DV',' 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)]; 3047830 27 - la détermination de l'angle /3 01 comme étant l'angle le' correspondant au maximum maxi., Mn,'des Mn,' calculés. Le procédé calcule alors en 504-6 la vitesse moyenne Vn des nuages, en pixels par 5 seconde, en se fondant sur la relation (2). La portion supérieure de la figure 21 illustre graphiquement la rotation de la matrice DV, illustrée à la figure 20 pour un angle /3v. La figure 22 illustre l'évolution du maximum Mn,' en fonction de l'angle len. On constate que le maximum de la courbe mi,,'(/3o 10 présente un maximum marqué correspondant à l'angle &opt,. La figure 23 illustre la détermination manuelle par un opérateur de l'angle /3 opta,. qui trace les pentes des principales structures de la figure 20. Le rapporteur à droite de la figure 23 montre que l'angle flopt,'obtenu par le procédé selon l'invention est cohérent avec une mesure de celui-ci faite directement sur l'image des distributions DVkopt,'. 15 En outre, l'écart de l'angle /3 opta,. 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 Dkoptan. 20 C.S. 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. 25 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 30 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 35 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). 3047830 28 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 DVkopta, 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 5 DVkopt,'), il 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 DVkopta, et de la vitesse moyenne Vn de déplacement des nuages. 10 Par exemple, l'estimation de l'ensoleillement consiste à détecter, en 506-2, les valeurs significativement non nulles de la dernière distribution DVkopta, 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 15 l'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 20 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 25 à 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 T,, en décalant temporellement les figures 25A-G de manière correspondante (de 50 / 5,67 = 8,81 pas de temps pour la Figure 25A, de 100 / 5,67 = 17,6 pas de temps pour la Figure 25B, de 150 / 30 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 35 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 3047830 29 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 5 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 10 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 Of , 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 15 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 calibration antérieure, il est déterminé que l'appareil a bougé si la différence entre les 20 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 {/A/7" est 25 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 chaque image IN': est calculée en fonction de l'image elle-même. Notamment, le calcul 30 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 comprend tout d'abord le calcul de la matrice de luminosité aci, et de la matrice 35 = IR(INnc) - B(INnc)1. 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 3047830 30 nvnc,R-B et monm balayage pixel par pixel des matrices M'in /Lcn est ensuite 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 5 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 le pixel est plus gris que bleu (i.e. illnc'R-B (a, b) REVENDICATIONS1. Procédé d'estimation de la direction moyenne de déplacements d'objets en mouvement dans un flux d'images numériques {MO d'une scène comprenant du ciel, les images numériques dudit flux étant munies d'un axe de référence Are f, le procédé comprenant : - l'acquisition du flux d'images numériques {MO ; - pour chaque image IN?, du flux d'images numériques {MO o le calcul d'une image segmentée INe, identifiant les objets en mouvement dans l'image INn; 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 A' f o pour chaque zone Z km, dudit ensemble de zones {Zka,} , le calcul d'une distribution DDk,i, des valeurs des pixels de ladite image segmentée INe, 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 km; - pour chaque direction d'analyse Dk, le calcul d'un critère de variabilité CVka, de distributions {DDka,_p, DDka,_p+i, , DDk,72_1, DDka,} calculées pour des zones Z km, associées à ladite direction Dk; et - la détermination de la direction moyenne Dkoptal 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é CV/c,' est optimum. 2. Procédé selon la revendication 1, caractérisé en ce que l'acquisition du flux d'image numériques {/Nn} 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 DDk,', le redressement des images les images numériques {/Nn} . 3. Procédé selon la revendication 1 ou 2, caractérisé en ce que 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 {/Nn} , et en 3047830 34 5 4. 10 5. 15 6. 20 7. 25 30 ce que le procédé comprend le calcul de l'ensemble de zones d'image {Zka,} pour chaque image I A In du flux d'images numériques {MO , 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 rectangles {Rka,} avec l'image I Ni,. Procédé selon la revendication 3, 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. Procédé selon la revendication 4, caractérisé en ce le petit côté de chaque rectangle Rical de l'ensemble de rectangles {Rk,72} est au moins égal au diamètre du halo Di,' du disque solaire. Procédé selon la revendication 4 ou 5, caractérisé en ce que 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. Procédé selon la revendication 6, caractérisé en ce que le modèle astronomique est réglé pour une position prédéterminée (Of , 0 , 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 en ce que le procédé comprend 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 {/N79 comprenant chacune le disque solaire ; - pour chaque image In dudit flux {/N79 , le calcul d'une position Prf du centre du disque solaire dans l'image I n en analysant l'image In ; - le calcul de la position du dispositif d'imagerie (Of , 0 , Azr) comme solution d'un problème d'optimisation selon a relation : /Pi 1 (Pf (Bi, 02, AZ)) (0f, 0 , Azr) = ar gmin (01,02,Az)612 vsnc:Pn (Bi, 02, AZ) 35 où D. est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart entre deux vecteurs de 3 3047830 35 8. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que le critère de de variabilité CVka, pour la direction Dk est calculé en fonction des distributions IDDk,n_p, DDk,n_p+i, , DDk,n1 associées aux p > 2 dernières images numériques acquises. 9. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que le critère de variabilité CVk,' pour la direction Dk est calculé selon la relation : CVkal = Dpk,j11 j=n-p+1 5 10 15 10. 20 11. 25 30 12. 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 (n-p)lème image dans le flux d'images numériques {/Nn} , avec p strictement inférieur à n, k référençant la klème direction Dk, et 11 V II est une norme du vecteur V de l'espace l'Zi , I étant la dimension des distributions DDk,n. Procédé selon l'une quelconque des revendications, 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,n, ledit masque est appliqué aux images du flux d'images numériques {/Nn} et/ou aux images segmentées {Min . 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 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 Dkopt,n. Procédé selon la revendication 11, 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. 3047830 36 13. Procédé selon la revendication 11 ou 12, caractérisé en ce que le calcul de ladite vitesse moyenne Vn comprend : pour chaque image segmentée INn, le calcul d'une distribution DVkopt,n des valeurs des pixels de l'image segmentée INe, appartenant à la zone Zkoptal 5 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 ; 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 10 flux jusqu'à la nième 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. 14. Procédé selon la revendication 13, caractérisé en ce que le calcul de l'inclinaison 15 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 DV, ; le calcul d'un histogramme H istnx pour chaque matrice DV,,' obtenue par rotation et le calcul du maximum Mn,' de cet histogramme ; 20 la détermination d'inclinaison comme étant égal à l'angle &opta,. correspondant au maximum des maxima Mn,' calculés. 15. Système d'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 25 ciel, les images numériques dudit flux étant munies d'un axe de référence Are f, le système comprenant : - une caméra pour l'acquisition du flux d'images numériques {INn} de ladite scène; une unité de traitement d'information configurée pour : 30 o pour chaque image INn du flux d'images numériques {INn} ^ le calcul d'une image segmentée INe, identifiant les objets en mouvement dans l'image INn; ^ le calcul d'un ensemble de zones d'image {Zka2} respectivement associées à des directions d'analyse {Dk} formant des angles {ak} 35 différents prédéterminés avec l'axe de référence A' f 3047830 37 ^ pour chaque zone Z kat dudit ensemble de zones {Zka,} , le calcul d'une distribution DDka, des valeurs des pixels de ladite image segmentée INn compris dans ladite zone Z kat, la distribution DDka, étant calculée le long d'une direction normale à la direction d'analyse 5 Dk de ladite zone Z km; o pour chaque direction d'analyse Dk, le calcul d'un critère de variabilité CVk,i, de distributions {DDka,_p, DDka,_p+i, , DDk,72_1, DDka,} calculées pour des zones Z kat associées à ladite direction Dk; et o la détermination de la direction moyenne Dkoptal de déplacement des 10 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. 16. Système selon la revendication 15 configuré pour mettre en oeuvre un procédé selon 15 l'une quelconque des revendications 2 à 14.