Brevet FR3047831A1 — publié le 20170818

PROCEDE DE DETERMINATION DE LA VITESSE 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 vitesse moyenne de déplacement d'objets dans la scène. L'invention trouve particulièrement application dans la détermination de la vitesse des nuages dans un ciel. La vitesse ainsi déterminée peut notamment être utilisée pour l'estimation de l'ensoleillement 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 3047831 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. 3047831 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 vitesse 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 vitesse moyenne 14, de déplacement d'objets en mouvement selon une direction Dkoptal dans un flux d'image 35 numérique {/Nn} d'une scène comprenant du ciel, comprenant : - l'acquisition du flux d'images numériques {/Nn} ; - l'estimation de la direction Dkoptal en fonction du flux d'images numériques {/Nn} - pour chaque image Mn du flux d'images numériques {/Nn} 3047831 4 o le calcul d'une image segmentée Me, identifiant les objets en mouvement dans l'image Mn; o le calcul d'une distribution DVkopt,n des valeurs des pixels d'au moins une zone Zkopt,n de l'image segmentée /Ne., la distribution DVkopt,n étant calculée le long 5 de la direction Dkopt,n 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,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; 10 et le calcul de la vitesse moyenne Vn en fonction de l'inclinaison calculée. En d'autres termes, les inventeurs ont constaté qu'il apparait naturellement des structures affines dans la matrice DVn = [DVkopt,n_qDVkopt,n_q+1 DVkopt,n_1DVkopt,n]. Ces structures affines et parallèles correspondent au déplacement moyen des objets dans le 15 flux d'image et leur pente correspond à la vitesse moyenne des objets selon la direction moyenne de déplacement. Le calcul de la vitesse moyenne est ainsi robuste et se limite à la détermination d'une pente. La direction de déplacement peut quant à elle être estimée par toute méthode appropriée connue de l'état de la technique, mais avantageusement selon l'estimation décrite ci-après. 20 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 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° ; 25 et préalablement au calcul des distributions DDk,n, le redressement des images les images numériques {/Nn} . 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 30 scène visible, et donc la portion visible du ciel. Selon un mode de réalisation, le calcul de l'inclinaison de la au moins une structure affine comprise dans une matrice DVn comprend : l'application de différentes rotations d'angles leu prédéterminés à la matrice DVn ; 35 le calcul d'un histogramme Histn,' pour chaque matrice DVn,n obtenue par rotation et le calcul du maximum Mn,' de cet histogramme ; la détermination d'inclinaison comme étant égal à l'angle lenopt,n correspondant au maximum des maxima Mn,' calculés. 3047831 5 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. Selon un mode de réalisation, l'estimation de la direction Dkoptal de déplacements des 5 objets en mouvement dans le flux d'images numériques {/Nn} d'une scène, comporte : pour chaque image Mn du flux d'images numériques {/Nn} 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 un axe de référence Are f des images numériques: 10 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 15 distributions {DDk,'_ DD k,n-p+1, --- D Dk,n_i, DDk,n} calculées pour des zones Zka, associées à ladite direction Dk; et la détermination de la direction Dkopt,n 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. 20 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 25 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 30 ressources ou une calibration, le procédé selon l'invention s'avère cependant très robuste à celle-ci de sorte qu'une segmentation peu précise des images est suffisante. Enfin, il a été constaté que le procédé est robuste également au processus optique de formation des images et renvoie des résultats très satisfaisants, même sur des images réalisées à partir d'un fisheye, y compris lorsqu'on leur applique un redressement sans traitement des 35 aberrations liées à ce type d'optique. 3047831 6 Selon un mode de réalisation, la scène est acquise à l'aide d'un dispositif d'imagerie comprenant un système optique apte à former une image de la scène sur un plan focal et un capteur numérique plan disposé sur le plan focal du système optique pour produire le flux d'images numérique {INn} , et le procédé comprend le calcul de l'ensemble de 5 zones d'image {Zka,} pour chaque image INn du flux d'images numériques {INn} , 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 10 - le calcul des zones {Zka,} comme étant égales respectivement à l'intersection des rectangles {Rka,} avec l'image INn. En particulier, la scène comprend du ciel, et le point prédéterminé du plan focal est la position du centre du disque solaire sur le plan focal du système optique. 15 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. 20 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 25 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 30 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. Plus particulièrement, le modèle astronomique est réglé pour une position prédéterminée 35 (61f, Oz, Azr) du dispositif d'imagerie, définie par la position (61f, 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 de la position prédéterminée (Of , Azr) du dispositif d'imagerie, comportant : 3047831 7 - 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 In en analysant l'image In ; 5 - le calcul de la position du dispositif d'imagerie (Of , Oz, Azr) comme solution d'un problème d'optimisation selon la relation : ( lpf.\ 131 (01, 02, Az) Vnc, / Pns. (01,02, Az) (Of , 02, Azr) = argmin (01,02,Az)612 où D. est un ensemble de recherche, et II II est une fonction de coût quantifiant un écart 10 entre deux vecteurs de II: 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, 15 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 = IR(INi2c) - 20 B(II\17,c) 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 Pnc dans la matrice I Lcn. En particulier, le calcul d'une position Prf du centre du disque solaire dans l'image In 25 comprend : le calcul de zones de pixels de la matrice /Ni2c'R-B candidats 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 30 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 particulier, le critère de de variabilité CVka, pour la direction Dk est calculé en fonction des distributions {DDk,'_p, DDk,72_p+1, DDk,72_1, DDk,72} 35 associées aux p > 2 dernières images numériques acquises. 3047831 8 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 5 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)iè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 II VII est une norme du vecteur V de l'espace 1 , I étant la dimension des distributions DDk,n. 10 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 15 distribution permet de détecter ce minimum avec précision. Selon un mode de réalisation, le procédé comporte le calcul d'un masque d'exclusion des éléments immobiles de la scène, et préalablement au calcul des distributions DDk,n, ledit masque est appliqué aux images du flux d'images numériques {/Nn} et/ou aux images 20 segmentées {/Nin . En d'autres termes, en excluant les objets de la scène qui ne comportent aucune information signifiante (arbres, reliefs, bâtiments, etc...), le traitement se concentre uniquement sur le ciel et voit donc sa précision renforcée. Selon un mode de réalisation, le procédé comprend en outre la prédiction de la position 25 d'un objet en mouvement en fonction de ladite vitesse Vn et de ladite direction Dkopt,n. 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. En particulier, le procédé comprend la 30 prévision de l'ensoleillement au-dessus d'un parc photovoltaïque, en particulier la chute ou l'accroissement d'ensoleillement. En d'autres termes, le procédé exploite avantageusement le fait que la direction de déplacement est déjà connue pour calculer la vitesse moyenne selon cette direction. 35 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 3047831 9 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. 5 L'invention a également pour objet un système d'estimation de la vitesse moyenne 142 de déplacement d'objets en mouvement selon une direction Dkoptal dans un flux d'images numériques {/Nn} d'une scène comprenant du ciel, le système comprenant : une caméra pour l'acquisition du flux d'images numériques {/Nn} de ladite scène; une unité de traitement d'information configurée pour : 10 o l'estimation d'une direction Dkoptal en fonction du flux d'images numériques {/Nn} ; o pour chaque image Mn du flux d'images numériques {/Nn} ^ le calcul d'une image segmentée Me, identifiant les objets en mouvement dans l'image Mn; 15 ^ le calcul d'une distribution DVkopt,n des valeurs des pixels d'au moins une zone Zkopt,n de l'image segmentée /Ne., la distribution DVkopt,n étant calculée le long de la direction Dkoptal déterminée ; o le calcul de l'inclinaison d'au moins une structure affine comprise dans une matrice DVn = [DVkopt,n_qDVkopt,'_q+1 DVkopt,n_iDVkopt,n] issue de la 20 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; o et le calcul de la vitesse moyenne Vn en fonction de l'inclinaison calculée. 25 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 vitesse moyenne Vn de déplacement d'objets 30 en mouvement selon une direction Dkoptal dans un flux d'image numérique {/Nn} d'une scène comprenant du ciel, le procédé comprenant : pour chaque image Mn du flux d'images numériques {/Nn} o le calcul d'une image segmentée Me, identifiant les objets en mouvement dans l'image Mn; 35 o le calcul d'une distribution DVkopt,n des valeurs des pixels d'au moins une zone Zkopt,n de l'image segmentée /Ne., la distribution DVkopt,n étant calculée le long de la direction Dkoptm déterminée ; 3047831 10 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)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; 5 et le calcul de la vitesse moyenne Vn en fonction de l'inclinaison calculée. En particulier, le support comprend des instructions pour la mise en oeuvre d'un procédé du type précité. 10 BREVE DESCRIPTION DES FIGURES L'invention sera mieux comprise à la lecture de la description qui va suivre, donnée uniquement à titre d'exemple, et faite en relation avec les dessins annexés dans lesquels : la figure 1 est une vue schématique d'un système de prédiction d'ensoleillement 15 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 ; la figure 3 est une image, transformée en noir et blanc, d'une image acquise par le 20 système de la figure 1 ; la figure 4 est une image issue d'une segmentation, e.g. de type Otsu, de la figure 3 pour la détermination d'objets non météorologiques ; la figure 5 est un profil illustrant la localisation d'une portion de ciel dans l'image de la figure 3, issu de la segmentation de la figure 4 ; 25 la figure 6 est une image numérique d'un masque d'exclusion obtenu à partir du profil de la figure 5 ; la figure 7 est une image numérique d'un masque d'exclusion obtenu par application d'un masque de dilation circulaire à l'image de la figure 6 ; la figure 8 est une image numérique illustrant l'application du masque de la figure 7 30 à l'image numérique de la figure 3 ; la figure 9 est une image numérique illustrant l'application d'un masque comprenant en outre le disque solaire à l'image numérique de la figure 3 ; la figure 10 est une image numérique illustrant une segmentation pour la détection des nuages dans l'image de la figure 9 ; 35 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 ; 3047831 11 - les figures 13A-E sont des images illustrant des flux de distributions de déplacement pour différents angles de détection ; - la figure 14 est une vue schématique illustrant graphiquement le calcul d'un flux de distribution ; 5 - les figures 15 à 17 sont des tracées illustrant différents critères de variabilité de flux en fonction de l'angle de détection ; - la figure 18 illustre une comparaison de la détection automatique de la direction moyenne de déplacement selon l'invention et une détection manuelle de cette direction ; 10 - la figure 19 est une vue schématique illustrant graphiquement le calcul d'une distribution de vitesse selon l'invention ; - la figure 20 est une image illustrant graphiquement des structures affines d'un flux de distributions de vitesse ; - la figure 21 est une image illustrant une rotation des structures des figures de la 15 figure 20 ; - la figure 22 est un tracé d'un critère d'une distribution de vitesse en fonction d'un angle de rotation de cette distribution; - la figure 23 illustre la comparaison entre la détermination automatique de l'angle de rotation correspondant à la vitesse de déplacement des nuages et la détermination 20 manuelle de cet angle ; - la figure 24 est une image illustrant une partie du flux de distribution de vitesse ; - les figures 25A-G sont des tracés illustrant des échéances temporelles d'arrivée de nuage sur le disque solaire pour différentes distances en pixel du soleil. - les figures 26 et 27 sont des tracées superposant les instants estimés d'arrivée des 25 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. 30 DESCRIPTION DETAILLEE DE L'INVENTION A. NOTATIONS Dans la suite du document, les conventions mathématiques classiques de l'algèbre 35 linéaire sont appliquées pour des matrices réelles de 11AxB où A et B sont des entiers positifs, avec notamment pour des matrices M et N deAxB M = (rna,b)aE[1,A],bE[1,B] 3047831 12 lm = (Ima,b1) aE[1,A],bE[1,B] 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 5 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 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 pixel de la matrice N en fonction du masque M peut être réalisée de nombreuses manières 10 connues en soit, la notion de masque étant classique en traitement d'images. Par convention, dans la demande la « matrice » M N est représentée graphiquement par l'image de N modifiée en réglant ses pixels en noir lorsqu'ils correspondent aux pixels du masque M de valeur 1. 15 B. MA 1 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 photographique ou une caméra 12, par exemple de type reflex, muni d'un capteur 20 numérique CMOS ou CCD 14 pour l'acquisition d'image dans le rayonnement visible, par exemple au format 35mm, et d'un objectif 16 à très grand angle de champ supérieur à 100° sur la diagonale du capteur 14, avantageusement de 180° (e.g. fisheye). L'objectif produit ainsi sur le capteur 14 une image optique circulaire de l'espace devant l'objectif 16 et donc une image comprenant l'ensemble du ciel visible depuis la position du 25 dispositif 10. Le capteur 14 est agencé sur le plan focal 18 de l'objectif 16 et centré sur l'axe optique 20 de ce dernier. L'appareil 12 comprend par ailleurs les circuits 22 nécessaires au pilotage de l'acquisition des images et le prétraitement des données issues du capteur (e.g. dématriçage) pour produire une image à partir des signaux générés depuis les photosites sensibles du capteur 14, usuellement un processeur « DSP » (pour « digital 30 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 35 la focale 8mm, piloté par un intervallomètre externe de référence DigiSnap 2700 de la société Harbonics. 3047831 13 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) 5 ou dans un format de lecture (e.g. jpeg), et pour leur appliquer un traitement mis en oeuvre par ordinateur selon l'invention. L'unité 24 est par exemple un ordinateur personnel, une tablette, un smartphone, ou plus généralement tout système à base de microprocesseur(s), notamment de type DSP 10 (« 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 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 15 mémorisation des résultats des calculs intermédiaires et finaux. 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 20 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 photovoltaïque est gérée. 25 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 de manière non régulière avec un pas d'acquisition non constant dans le temps. 30 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 composant les matrices sont codées sur un nombre bits prédéterminé, par exemple 8 bits, 35 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). 3047831 14 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, 5 l'invention s'applique à un nombre de bits quelconque. C. TRAI I EMENT INFORMATIQUE Il va à présent être décrit, en relation avec les organigrammes des figures 2A-2E, un 10 procédé selon un mode de réalisation de prévision de la production d'électricité d'un parc photovoltaïque, et plus particulièrement la chute d'ensoleillement sur le parc photovoltaïque 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. 15 C.1. Calibration Le procédé débute par une phase 40 de positionnement et de calibration du dispositif d'imagerie 10. 20 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, 25 la latitude, de l'heure et du jour de l'année. Ce modèle est par exemple le « solar position calculator » de la NOAA (« National Oceanic and Atmospheric Administration »); - un modèle de la formation d'image sur le capteur 14 par l'objectif 16. Pour le fisheye, il s'agit par exemple du modèle décrit dans le document « Geometric 30 modelling and calibration of fisheye lens camera system », Schwalbe E., In Proc. 2nd Panaramic Photogrammetry Worshop, Int. Archives of Photogrammetry and Remote Sensing (vol. 36, No. 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. 35 3047831 15 En connaissant la position (01, 02, 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 5 produites par l'appareil 12. Cette combinaison est ci-après désignée par les termes « modèle astronomique optique ». En 402, un opérateur positionne l'appareil 12 sur un site choisi avec l'objectif visant le ciel. La position de l'appareil est ainsi défini par trois paramètres (01, 02,Az), le couple 10 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 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 15 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 mémorisées dans l'unité 24. Evidemment, il est recherché à positionner le dispositif dans une zone dégagée de 20 manière à maximiser la portion de ciel dans les images. Toutefois, le dispositif peut être implanté dans une zone très contrainte en termes d'éléments non pertinents présents dans le champ de l'appareil (ou « objet non météorologiques », ou « onm », par exemple des arbres, des constructions, des pilonnes, etc...), l'invention permettant en effet une prévision robuste y compris sur des portions minimes de ciel et y compris avec le soleil 25 caché. 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 (0f, 0 , Azr) différente de la position théorique (01, 02, Az). La position 30 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 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 35 dans les images produites par l'appareil 12. 3047831 16 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 5 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. 10 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é 15 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 20 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 25 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 histograms », de Otsu N, Automatica, 11(285-296), 23-27, 1975. Ce type 30 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) = argmin (6(0, s1) + o-(s1 + 1, s2) + o- (s2 + 1, s3) + o- (s3 + 1,255)) (sl,s2,53)E [0,255] 3 35 3047831 17 où 6(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 3047831 24 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 3047831 25 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 3047831 26 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 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 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 MW,u(J3is) 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. 3047831 27 En outre, l'écart de l'angle &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 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 3047831 28 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 3047831 29 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 vitesse moyenne Vn de déplacement d'objets en mouvement selon une direction Dkopt,ii dans un flux d'images numériques {INn} d'une scène comprenant du ciel, comprenant : - l'acquisition du flux d'images numériques {INn} - l'estimation de la direction moyenne Dkopt,r1 en fonction du flux d'images numériques {INn} ; - pour chaque image INn du flux d'images numériques {INn} : 10 o le calcul d'une image segmentée In identifiant les objets -en mouvement dans l'image INn; le calcul d'une distribution DVkopt,n des valeurs des pixels d'au moins une zone Zkopt,n de l'image segmentée INL la distribution DVkopt,, étant calculée le long de la direction moyenne Dkopt,n déterminée ; 15 le calcul de l'inclinaison d'au moins une structure affine comprise dans une matrice DVn =[D17kopt,n-qpi/kopt,n-q+1.--DV_ k opt, _-_DV_k opt, _1 issue de la concaténation des distributions DVkopt,' des images entre la (n-q)1eme image du flux jusqu'à la nième du flux d'images numériques {INn} , où q est un entier supérieur ou égal à 1; 20 - et le calcul de la vitesse moyenne Vn en fonction de l'inclinaison calculée. 2. Procédé selon la revendication 1, caractérisé en ce que l'acquisition du flux d'image numériques {INn} comprend : l'acquisition de la scène par un dispositif de prise de vue ayant un angle de 25 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 - le redressement des images numériques {INn} 30 3. Procédé selon la revendication 1 ou 2, caractérisé en ce que le calcul de l'inclinaison de la au moins une structure affine comprise dans une matrice DVn comprend : - l'application de différentes rotations d'angles Ai prédéterminés à la matrice DVn; 35 le calcul d'un histogramme Histnx pour chaque matrice DVn,n obtenue par rotation et le calcul du maximum Mn u de cet histogramme ; - la détermination d'inclinaison comme étant égal à l'angle pito ,' correspondant au maximum des maxima Mn,' calculés. 3047831 33 4. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que l'estimation de la direction Dkopt,n de déplacements des objets en mouvement dans le flux d'images numériques {INn} d'une scène, comporte : - pour chaque image INn du flux d'images numériques t/N7,1 o le calcul d'un ensemble de zones d'image {Zk,,,} respectivement associées à des directions d'analyse {Dk} formant des angles {alc} différents prédéterminés avec un axe de référence Are! des images numériques: o pour chaque zone Zkii dudit ensemble de zones {Zk,'} , le calcul d'une distribution DDk,' des valeurs des pixels de ladite image segmentée in compris dans ladite zone Zk,n, la distribution DDk;' étant calculée le long d'une direction normale à la direction d'analyse Dk de ladite zone Zkn; pour chaque direction d'analyse Dk, le calcul d'un critère de variabilité CVk,' de distributions {DDk,n_p, DDk,n_p+i, , DDk,'} calculées pour des zones Zkn associées à ladite direction Dk; et la détermination de la direction Dkopt,n de déplacement des objets dans le flux d'images numériques {INn} comme étant égale à la direction d'analyse pour laquelle le critère de variabilité CVk,,, est optimum. 5. Procédé selon la revendication 4, 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 {INn} , et en ce que le procédé comprend le calcul de l'ensemble de zones d'image {Zk,,,} pour chaque image INn du flux d'images numériques {INn} , ledit calcul comprenant : le calcul d'un ensemble de rectangles {Rk,n} 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 14,,J comme étant égales respectivement à l'intersection des rectangles {Rk,n} avec l'image INn. 6. Procédé selon la revendication 5, 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. 3047831 34 7. Procédé selon la revendication 6, caractérisé en ce le petit côté de chaque rectangle Rkm de l'ensemble de rectangles {Rk } est au moins égal au diamètre du halo De du disque solaire. 8. Procédé selon la revendication 6 ou 7, 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. 10 9. Procédé selon la revendication 8, caractérisé en ce que le modèle astronomique est réglé pour une position prédéterminée (Of ,6,I,Azr) du dispositif d'imagerie, définie par la position (Of ,01) 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 15 prédéterminée (61f, 0Z, Azr) du dispositif d'imagerie, comportant : l'acquisition d'un flux d'images numériques de calibration horodatées (IN79 comprenant chacune le disque solaire ; pour chaque image INn dudit flux {//g) , le calcul d'une position Pif du centre du disque solaire dans l'image In en analysant l'image In ; 20 le calcul de la position du dispositif d'imagerie (01, 01, Azr) comme solution d'un problème d'optimisation selon la relation : (Of , Azr) = ar gmin (Pic. (Pi' (611, 02, Az) (91,e2,Az)En Prf Pr'1.5 (01,. 612, Az) ._ où D. est un ensemble de recherche, et II II est une fonction de coût quantifiant 25 un écart entre deux vecteurs de 1[83. 10. Procédé selon l'une quelconque des revendications 4 à 9, caractérisé en ce que le critère de de variabilité CVk,' pour la direction Dk est calculé en fonction des distributions {DDk,n_p, DDk,n_p÷ 1, , DDk,n} associées aux p 2 30 dernières images numériques acquises. 11. Procédé selon l'une quelconque des revendications 4 à 10, caractérisé en ce que le critère de variabilité CVk,n pour la direction Dk est calculé selon la relation CVkm, = Il DDk - DDki II j--=h-p+1 3047831 35 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 dème image et la (n-p)1ème image dans le flux d'images numériques, {INn} , avec p strictement inférieur à n, k référençant la kième direction Dk, et 11 V II est une norme du vecteur V de l'espace Ri , I étant la dimension des distributions DDk,n. 12. Procédé selon l'une quelconque des revendications 4 à 11, 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 10 appliqué aux images du flux d'images numériques {INn} et/ou aux images segmentées {Min . 13. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il comprend en outre la prédiction de la position d'un objet en mouvement en 15 fonction de ladite vitesse moyenne lin et de ladite direction Dkoptm. 14. Procédé selon la revendication 13, 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. 20 15. Système d'estimation de la vitesse moyenne Vn de déplacement d'objets en mouvement selon une direction Dkopt,n dans un flux d'images numériques {INn} d'une scène comprenant du ciel, le système comprenant : une caméra pour l'acquisition du flux d'images numériques {INn} de ladite 25 scène; - une unité de traitement d'information configurée pour : o l'estimation d'une direction Dkopt,n en fonction du flux d'images numériques {INn} o pour chaque image INn du flux d'images numériques {Mn} : 30 ^ le calcul d'une image segmentée In identifiant les objets en mouvement dans l'image INn; ^ le calcul d'une distribution Dlikopt,n des valeurs des pixels d'au moins une zone Zkopt,n de l'image segmentée In, la distribution DVkopt,n étant calculée le long de la direction Dkoptm déterminée ; 35 3047831 36 o le calcul de l'inclinaison d'au moins une structure affine comprise dans une matrice DV' = [DVkopt,n_eVkopt,'_q+1 DVkopt,,,_1DVIcopt,,,] issue de la concaténation des distributions DVkopt,' des images entre la (n-q)ième image du flux jusqu'à la nié' du flux d'images numériques (Mn} , où q est un entier supérieur ou égal à 1; o et le calcul de la vitesse moyenne 1Vn en fonction de l'inclinaison calculée. 16. Système selon la revendication 15 configuré pour mettre en oeuvre un procédé selon l'une qpelconque des revendications 2 à 14.