FR3047831A1 - Procede de determination de la vitesse de deplacement d'objets dans une scene - Google Patents

Abstract

Un procédé d'estimation de la vitesse moyenne Vn de déplacement d'objets en mouvement dans un flux d'image numérique {INn} d'une scène comprenant du ciel, comprend : - l'acquisition du flux {INn} ; - l'estimation d'une direction moyenne Dkopt,n de déplacement des objets; pour chaque image INn du flux : ○ le calcul d'une image segmentée INn 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 INSn, la distribution DVkopt,n étant calculée le long de la direction moyenne 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] concaténant des distributions DVkopt,n entre la (n-q)ième image et la nième image du flux, 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.

Description

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 +lh à +48h. L’advection d’images satellite permet notamment de prévoir la nébulosité aux échéances les plus courtes (jusqu’à +5h), tandis que l’adaptation statistique des modèles météorologiques permet de calculer des prévisions aux échéances plus lointaines.

Cependant, à l’échelle infra-journalière, la production d’énergie sur les parcs photovoltaïques peut chuter brutalement en quelques dizaines de secondes lors du passage de nuages au-dessus du site d’implantation. La reprise est tout aussi soudaine. Cette intermittence est problématique quand l’énergie du parc est utilisée en mode stockage/déstockage en association avec des batteries ou d’autres moyens de production d’énergie renouvelable (centrales hydroélectriques). Les mécanismes d’asservissement sur ces temps de réactions relativement courts peuvent cependant être optimisés si la production d’énergie par les panneaux photovoltaïques est anticipée de quelques minutes à quelques heures. En outre, le couvert nuageux a non seulement une influence sur l’amplitude mais également sur la composition du rayonnement global. La présence de nuages modifie sensiblement l’importance relative des rayonnements diffus et direct. Alors que par ciel clair le rayonnement diffus représente aux mieux 10 à 15% du rayonnement global, les nuages constituent un facteur de diffusion supplémentaire. Par ciel couvert de nuages élevés, l’essentiel du rayonnement global est attribuable au rayonnement diffus.

Or, les techniques à base de modèle météorologique et d’images satellite ne permettant de traiter ni les phénomènes locaux, ni les phénomènes à courte échelle de temps, ni de tenir compte du couvert nuageux. C’est pourquoi l’approche le plus souvent adoptée pour appréhender et anticiper ces phénomènes au voisinage immédiat d’un parc photovoltaïque consiste à acquérir des images du ciel au-dessus du parc depuis le sol, à les traiter pour déterminer la course des nuages dans le ciel et tenter de prévoir sur cette base l’ensoleillement au-dessus du parc photovoltaïque.

Toutefois, cette approche pose des problèmes tant au niveau du matériel mis en œuvre que du traitement d’images réalisé sur le flux d’images acquises. Par exemple acquérir l’image du ciel en une seule prise nécessite l’emploi d’une optique ayant un très grand angle de champ, par exemple de type « fish-eye » qui permet d’obtenir des images hémisphériques englobant la totalité du ciel visible depuis la position de prise de vue. Cependant, la capacité à prendre des photos de très grand angle s’accompagne de nombreux défauts optiques : des distorsions géométriques importantes, des aberrations chromatiques ainsi que d’autres défauts optiques qui mettent ultérieurement en difficulté le traitement d’image et/ou faussent son résultat. Outre la distorsion plus difficile à corriger que pour des optiques classiques, l’un des inconvénients majeurs de ce type d’optique est que la résolution de l’image est très mauvaise lorsque que l’on se rapproche de l’horizon. C’est pourquoi des optiques sophistiquées ont été conçue (e.g. des systèmes catadioptriques panoramique à base de miroir hémisphérique, des miroirs convexes à bande tournante, etc.). Outre leur complexité, ce type d’optique impose l’emploi de pièces mécaniques mobiles (e.g. d’un traqueur pour suivre la course du soleil ou d’un moteur pour entraîner en rotation des éléments optiques), ce qui renchérit le coût global du dispositif, et pose de plus des problèmes de fiabilité et de robustesse.

Concernant le traitement d’image, les techniques les plus répandues pour suivre (« tracking ») des objets en mouvement dans une scène sont de type « Particle Image Velocity » (« PIV »). Cette technique utilise une recherche de ressemblance entre deux images espacées temporellement d’un court laps de temps à partir desquelles il est extrait un champ de vecteurs locaux de déplacement instantané (i.e : direction de mouvement et amplitude de la vitesse). Si cette technique présente des avantages, notamment la capacité à distinguer les zones uniformes des zones fortement variables, elle présente cependant des désavantages nombreux. Notamment, une partie significative de l’image est très souvent non renseignée sur les vecteurs de déplacements et/ou est faussement identifiée comme exempte d’objet en déplacement. Ceci nécessite une calibration très fine du traitement et/ou une post-analyse complexe et coûteuse en paramètres seuils, et donc particulièrement complexe à calibrer. De plus une identification correcte des déplacements par cette technique suppose des déformations faibles des objets (nuages) d’une image à l’autre, ce qui impose des prises d’images à intervalle de temps rapproché. Or, les nuages sont par nature de forme très changeante, notamment parce qu’ils se déforment, parce que de nouveaux nuages apparaissent ou disparaissent, etc... De plus quand les nuages approchent du disque solaire (cas particulièrement intéressant puisqu’en interceptant le disque solaire ces nuages affectent la radiation solaire perçue au sol), leurs propriétés colorimétriques sont souvent très fortement modifiées. Ils deviennent notamment plus clairs dans l’image et mettent donc en difficulté la recherche de ressemblance. Au final, ce type de technique est donc peu robuste pour la détermination de la course des nuages et nécessite d’importantes ressources informatiques pour la mise en œuvre du traitement.

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 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 Vn de déplacement d’objets en mouvement selon une direction Dkoptn dans un flux d’image numérique {INn} d’une scène comprenant du ciel, comprenant : l’acquisition du flux d’images numériques {INn} ; l’estimation de la direction DkoptTl en fonction du flux d’images numériques {INn} ; pour chaque image INn du flux d’images numériques {INn} o le calcul d’une image segmentée IN% identifiant les objets en mouvement dans l’image INn, o le calcul d’une distribution DVk n des valeurs des pixels d’au moins une zone Zkopt,n de l’image segmentée INP, la distribution DVk n étant calculée le long de la direction Dk n déterminée ; le calcul de l’inclinaison d’au moins une structure affine comprise dans une matrice DVn = [DVkoptin_qDVkoptin-q+1... DV^^DV^n] issue de la concaténation des distributions DVk n des images entre la (n-q)ieme image du flux jusqu’à la nieme du flux d’images numériques {INn} , où q est un entier supérieur ou égal à 1; et le calcul de la vitesse moyenne Vn en fonction de l’inclinaison calculée.

En d’autres termes, les inventeurs ont constaté qu’il apparaît naturellement des structures affines dans la matrice DVn = [DVkoptin.qDVkoptin-q+1 ...DV^^DV^n]. Ces structures affines et parallèles correspondent au déplacement moyen des objets dans le flux d’image et leur pente correspond à la vitesse moyenne des objets selon la direction moyenne de déplacement. Le calcul de la vitesse moyenne est ainsi robuste et se limite à la détermination d’une pente. 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.

Selon un mode de réalisation, l’acquisition du flux d’image numériques {/„} comprend : - l’acquisition de la scène par une dispositif de prise de vue ayant un angle de champ supérieur ou égal à 100° selon la diagonale d’un capteur au format 35mm ou en équivalent au format 35mm, en particulier un angle de champ égal à 180° ; - et préalablement au calcul des distributions DDkn, le redressement des images les images numériques {INn} L’invention exploite ainsi avantageusement les optiques à très angle de champ, par exemple les fisheyes qui permettent en une seule prise de vue d’acquérir l’ensemble de la scène visible, et donc la portion visible du ciel.

Selon un mode de réalisation, le calcul de l’inclinaison de la au moins une structure affine comprise dans une matrice DVn comprend : - l’application de différentes rotations d’angles βη prédéterminés à la matrice DVn ; - le calcul d’un histogramme Histnu pour chaque matrice DVnu obtenue par rotation et le calcul du maximum Mnu de cet histogramme ; - la détermination d’inclinaison comme étant égal à l’angle βηορί,η correspondant au maximum des maxima Mnu calculés.

Ce calcul permet de déterminer de manière précise l’inclinaison à l’aide de calcul algébriques simples, et donc peu consommateurs de ressources informatiques.

Selon un mode de réalisation, l’estimation de la direction Dkoptiîl 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 {INn} o le calcul d’un ensemble de zones d’image {Zkn} 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 Aref des images numériques: o pour chaque zone Zkn dudit ensemble de zones {Zkn} , le calcul d’une distribution DDk n des valeurs des pixels de ladite image segmentée IN£ compris dans ladite zone Zk n, la distribution DDk n étant calculée le long d’une direction normale à la direction d’analyse Dk de ladite zone Zfcn; pour chaque direction d’analyse Dk, le calcul d’un critère de variabilité CVkn de distributions {DDkn_p, DDkn_p+1, ....DD^^.DDun} calculées pour des zones Zkn associées à ladite direction Dk, et la détermination de la direction DkoptU de déplacement des objets dans le flux d’images numériques {INn} comme étant égale à la direction d’analyse pour laquelle le critère de variabilité CVk n est optimum.

En d’autres termes, les inventeurs ont découvert que la zone orientée qui présente la variation temporelle minimale de sa distribution est celle dont la direction correspond à la direction moyenne de déplacement des objets dans une scène. Ce procédé est ainsi bien adapté à la course des nuages dans le ciel. 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 œuvre sont basiques (calcul de distributions et calculs des variations de distribution), et ne nécessitent donc pas de grandes ressources informatiques. De plus, si la segmentation peut nécessiter des ressources ou une calibration, le procédé selon l’invention s’avère cependant très robuste à celle-ci de sorte qu’une segmentation peu précise des images est suffisante. Enfin, il a été constaté que le procédé est robuste également au processus optique de formation des images et renvoie des résultats très satisfaisants, même sur des images réalisées à partir d’un fisheye, y compris lorsqu’on leur applique un redressement sans traitement des aberrations liées à ce type d’optique.

Selon un mode de réalisation, la scène est acquise à l’aide d’un dispositif d’imagerie comprenant un système optique apte à former une image de la scène sur un plan focal et un capteur numérique plan disposé sur le plan focal du système optique pour produire le flux d’images numérique {INn} , et le procédé comprend le calcul de l’ensemble de zones d’image {Zfcn} pour chaque image INn du flux d’images numériques {INn} , ledit calcul comprenant : le calcul d’un ensemble de rectangles {Æfcn} ayant chacun un petit côté centré sur un même point prédéterminé du plan focal, et ayant leurs grands côtés respectivement parallèles aux directions d’analyse {Dk} ; et le calcul des zones {Zfcn} comme étant égales respectivement à l’intersection des rectangles {Æfcn} avec l’image INn.

En particulier, la scène comprend du ciel, et le point prédéterminé du plan focal est la position du centre du disque solaire sur le plan focal du système optique.

En d’autres termes, le procédé permet de calculer directement le déplacement moyen dans un référentiel d’intérêt, par exemple centré sur l’image du soleil. Il est ainsi déterminé directement les nuages qui se déplacent vers le soleil et qui sont donc susceptible de masquer ultérieurement ce dernier.

Avantageusement, le petit côté de chaque rectangle Rkn de l’ensemble de rectangles {Æ/c,n} est au moins égal au diamètre du halo du disque solaire. Les images sont ainsi découpées en zones qui tiennent compte directement d’un résultat recherché ultérieurement, notamment l’occultation du disque solaire et de son halo par des nuages et donc la variation de l’ensoleillement.

Selon une variante avantageuse, la position du centre du disque solaire est calculée en fonction d’un modèle astronomique de l’azimut et de l’élévation du soleil dans le ciel définissant la position dudit centre sur le plan focal. En d’autres termes, il est possible de déterminer le déplacement des nuages vers le soleil y compris lorsque ce dernier est absent de l’image, par exemple lors qu’il est caché par un nuage, par un arbre, voire même complètement en dehors de l’image.

Plus particulièrement, le modèle astronomique est réglé pour une position prédéterminée (θ{, 02,Azr) du dispositif d’imagerie, définie par la position 0J) 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 (θ[, ,Azr) du dispositif d’imagerie, comportant : l’acquisition d’un flux d’images numériques de calibration horodatées {ΙΝξ} comprenant chacune le disque solaire ; pour chaque image IN£ dudit flux {ΙΝξ} , le calcul d’une position P£ du centre du disque solaire dans l’image IN£ en analysant l’image IN£ ; le calcul de la position du dispositif d’imagerie (θ[, ,rizr) comme solution d’un problème d’optimisation selon la relation :

où Ω est un ensemble de recherche, et || || est une fonction de coût quantifiant un écart entre deux vecteurs de M3.

En 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, automatiquement et avec précision cette erreur de positionnement.

Notamment, les images sont en couleur, et en ce que le calcul d’une position P£ du centre du disque solaire dans l’image IN£ comprend : le calcul d’une matrice de luminosité ILcn d’une matrice IN^R~B = \R(1N%) — B(IN%) | égale à la valeur absolue de la différence entre des canaux bleu et rouge de l’image IN£ ; et le calcul de la position P£ dans la matrice ILcn.

En particulier, le calcul d’une position P£ du centre du disque solaire dans l’image IN£ comprend : le calcul de zones de pixels de la matrice IN^R~Bcandidats au disque solaire, un pixel (a, h) étant jugé candidat si sa valeur Z£(a, b) dans la matrice ILcn est supérieure à un premier seuil prédéterminé et si sa valeur inn' (a, b) dans la matrice IN^R~B est inférieure à un second seuil prédéterminé ; et la sélection d’une zone candidate comme appartenant au disque solaire en fonction de sa proximité à une position prédéterminée.

Selon un mode de réalisation particulier, le critère de de variabilité CVkn pour la direction Dk est calculé en fonction des distributions {DDkn_p,DDkn_p+1, ...,DDkn_1, DDkn) associées aux p > 2 dernières images numériques acquises.

Selon un mode de réalisation, le critère de variabilité CVkn pour la direction Dk est calculé selon la relation :

expression dans laquelle j référence la position d’une image numérique dans le flux d’images numériques {INn} , n et n-p référence respectivement la nieme image et la (n-p)ieme image dans le flux d’images numériques {/7Vn} , avecp strictement inférieur àn,k référençant la kieme direction Dk, et ||E|| est une norme du vecteur V de l’espace R7 , I étant la dimension des distributions DDkn.

Les inventeurs ont en effet noté que le critère CVk n présente un minimum global en fonction de l’indice k et que ce minimum global correspond à la direction moyenne de déplacement. Le critère peut être plus au moins raffiné pour mettre en valeur ce minimum global mais même la simple somme de somme des valeurs de différences deux à deux de distribution permet de détecter ce minimum avec précision.

Selon un mode de réalisation, le procédé comporte le calcul d’un masque d’exclusion des éléments immobiles de la scène, et préalablement au calcul des distributions DDkn, ledit masque est appliqué aux images du flux d’images numériques {INn} et/ou aux images segmentées {IN%} . 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 d’un objet en mouvement en fonction de ladite vitesse Vn et de ladite direction Dkovt 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 {/JVn} ayant ladite vitesse et ladite direction. En particulier, le procédé comprend la prévision de l’ensoleillement au-dessus d’un parc photovoltaïque, en particulier la chute ou l’accroissement d’ensoleillement.

En d’autres termes, le procédé exploite avantageusement le fait que la direction de déplacement est déjà connue pour calculer la vitesse moyenne selon cette direction. Disposant d’une information importante, le calcul de la vitesse est donc facilité, tout type de calcul peut être utilisé à cette fin. Connaissant la direction et la vitesse moyenne, par exemple des nuages dans le ciel, il devient alors possible de faire des prévisions infra-journalières de l’ensoleillement, et donc de la production d’électricité par un parc photovoltaïque. L’invention a également pour objet un système d’estimation de la vitesse moyenne Vn de déplacement d’objets en mouvement selon une direction Dk 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 scène; - une unité de traitement d’information configurée pour : o l’estimation d’une direction DkoptTl en fonction du flux d’images numériques {/AU ; o pour chaque image INn du flux d’images numériques {INn} le calcul d’une image segmentée IN£ identifiant les objets en mouvement dans l’image INn; le calcul d’une distribution DVk n des valeurs des pixels d’au moins une zone ZkoptJl de l’image segmentée IN£, la distribution DVk n étant calculée le long de la direction Dk n déterminée ; o le calcul de l’inclinaison d’au moins une structure affine comprise dans une matrice DVn — ^DVkoptn_qDVkoptn_q+i... //ffe0pt,n-i//^kopt,nj issue de la concaténation des distributions DVk n des images entre la (n-q)ieme image du flux jusqu’à la nlème du flux d’images numériques {INn} , 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.

En particulier, le système est configuré pour mettre en œuvre un procédé du type précité. L’invention a également pour objet un produit programme d’ordinateur enregistré sur support informatique utilisable par un ordinateur comprenant des instructions pour l’exécution d’un procédé d’estimation de la vitesse moyenne Vn de déplacement d’objets en mouvement selon une direction Dk n dans un flux d’image numérique {INn} d’une scène comprenant du ciel, le procédé comprenant : - pour chaque image INn du flux d’images numériques {INn} o le calcul d’une image segmentée IN% identifiant les objets en mouvement dans l’image INn; o le calcul d’une distribution DVk n des valeurs des pixels d’au moins une zone Zkoptin de l’image segmentée IN£, la distribution DVk n étant calculée le long de la direction Dk n déterminée ; le calcul de rinclinaison d’au moins une structure affine comprise dans une matrice DVn = [DVkopt)n_qDVkopt)n_q+1 ... Wfcopt,„-iWfcopt,n] issue de la concaténation des distributions DVk n des images entre la (n-q)ieme image du flux jusqu’à la nieme du flux d’images numériques {INn} , où q est un entier supérieur ou égal à 1; et le calcul de la vitesse moyenne Vn en fonction de l’inclinaison calculée.

En particulier, le support comprend des instructions pour la mise en œuvre d’un procédé du type précité.

BREVE DESCRIPTION DES FIGURES L’invention sera mieux comprise à la lecture de la description qui va suivre, donnée uniquement à titre d’exemple, et faite en relation avec les dessins annexés dans lesquels : la figure 1 est une vue schématique d’un système de prédiction d’ensoleillement selon l’invention ; les figures 2A-2D sont des organigrammes d’un procédé de prévision d’ensoleillement prédiction de la production photovoltaïque mis en œuvre, par exemple mis en œuvre par le système de la figure 1 ; la figure 3 est une image, transformée en noir et blanc, d’une image acquise par le système de la figure 1 ; - la figure 4 est une image issue d’une segmentation, e.g. de type Otsu, de la figure 3 pour la détermination d’objets non météorologiques ; - la figure 5 est un profil illustrant la localisation d’une portion de ciel dans l’image de la figure 3, issu de la segmentation de la figure 4 ; - la figure 6 est une image numérique d’un masque d’exclusion obtenu à partir du profil de la figure 5 ; - la figure 7 est une image numérique d’un masque d’exclusion obtenu par application d’un masque de dilation circulaire à l’image de la figure 6 ; - la figure 8 est une image numérique illustrant l’application du masque de la figure 7 à l’image numérique de la figure 3 ; - la figure 9 est une image numérique illustrant l’application d’un masque comprenant en outre le disque solaire à l’image numérique de la figure 3 ; la figure 10 est une image numérique illustrant une segmentation pour la détection des nuages dans l’image de la figure 9 ; la figure 11 est une image illustrant un rectangle de détection selon l’invention ; - les figures 12A-C sont des figures illustrant différent type de zones de détection selon l’invention ; les figures 13A-E sont des images illustrant des flux de distributions de déplacement pour différents angles de détection ; - la figure 14 est une vue schématique illustrant graphiquement le calcul d’un flux de distribution ; - les figures 15 à 17 sont des tracées illustrant différents critères de variabilité de flux en fonction de l’angle de détection ; - la figure 18 illustre une comparaison de la détection automatique de la direction moyenne de déplacement selon l’invention et une détection manuelle de cette direction ; - la figure 19 est une vue schématique illustrant graphiquement le calcul d’une distribution de vitesse selon l’invention ; - la figure 20 est une image illustrant graphiquement des structures affines d’un flux de distributions de vitesse ; - la figure 21 est une image illustrant une rotation des structures des figures de la figure 20 ; - la figure 22 est un tracé d’un critère d’une distribution de vitesse en fonction d’un angle de rotation de cette distribution; - la figure 23 illustre la comparaison entre la détermination automatique de l’angle de rotation correspondant à la vitesse de déplacement des nuages et la détermination manuelle de cet angle ; la figure 24 est une image illustrant une partie du flux de distribution de vitesse ; les figures 25A-G sont des tracés illustrant des échéances temporelles d’arrivée de nuage sur le disque solaire pour différentes distances en pixel du soleil, les figures 26 et 27 sont des tracées superposant les instants estimés d’arrivée des nuages sur le soleil, la figure 27 illustrant en outre des mesures d’éclairement par un pyranomètre ; et les figures 28 et 29 sont des images illustrant le procédé de calibration de la position du dispositif de prise de vue selon l’invention.

Description detaillee de l’invention A. Notations

Dans la suite du document, les conventions mathématiques classiques de l’algèbre linéaire sont appliquées pour des matrices réelles de RAxB, où A et B sont des entiers positifs, avec notamment pour des matrices M et N de RAxB : M - (ma,ù)ae[ly4]jùe[1B] \M\ (\ma.b\)ae[l,A],be[l,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 ® IV signifie que seuls les pixels de la matrice N correspondant à des pixels du masque M ayant une première valeur sont pris en compte pour le traitement. Par exemple ma b = 1 si le pixel est pris en compte ou 0 s’il est exclu du traitement. De même, lorsqu’il est fait référence à l’objet M 0 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 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 TV modifiée en réglant ses pixels en noir lorsqu’ils correspondent aux pixels du masque M de valeur 1. B. MATERIEL·

En se référant à la figure 1, un dispositif 10 d’imagerie et de traitement d’image selon l’invention, directement implanté sur un parc photovoltaïque, comporte un appareil photographique ou une caméra 12, par exemple de type reflex, muni d’un capteur numérique CMOS ou CCD 14 pour l’acquisition d’image dans le rayonnement visible, par exemple au format 35mm, et d’un objectif 16 à très grand angle de champ supérieur à 100° sur la diagonale du capteur 14, avantageusement de 180° (e.g. fisheye). L’objectif produit ainsi sur le capteur 14 une image optique circulaire de l’espace devant l’objectif 16 et donc une image comprenant l’ensemble du ciel visible depuis la position du dispositif 10. Le capteur 14 est agencé sur le plan focal 18 de l’objectif 16 et centré sur l’axe optique 20 de ce dernier. L’appareil 12 comprend par ailleurs les circuits 22 nécessaires au pilotage de l’acquisition des images et le prétraitement des données issues du capteur (e.g. dématriçage) pour produire une image à partir des signaux générés depuis les photosites sensibles du capteur 14, usuellement un processeur « DSP » (pour « digital signal processor», ainsi qu’une horloge permettant d’horodater les images acquises. L’acquisition des images est réalisée avec une période d’acquisition Te, l’appareil comportant à cet effet un intervallomètre interne ou le dispositif comportant un intervallomètre externe d’une manière connu en soi. Par exemple, l’appareil 12 est un appareil Canon de référence 6D équipé d’un objectif Canon EF 8-15 mm f/4 L réglé sur la focale 8mm, piloté par un intervallomètre externe de référence DigiSnap 2700 de la société Harbonics.

Le dispositif 10 comporte par ailleurs une unité de traitement d’informations 24 connectée à l’appareil 12, par un câble ou par une communication sans fil, pour recevoir les images acquises par l’appareil 12, par exemple au format brut (après le dématriçage mais avant tout traitement informatique, correspondant au format « RAW » chez Canon) ou dans un format de lecture (e.g. jpeg), et pour leur appliquer un traitement mis en œuvre par ordinateur selon l’invention. L’unité 24 est par exemple un ordinateur personnel, une tablette, un smartphone, ou plus généralement tout système à base de microprocesseur s), notamment de type DSP (« digital signal processor »), à base de circuits de type FPGA, à base de circuits mixant ces types de technologie, etc., apte à mettre en œuvre l’invention. L’unité 24 est notamment pourvue de l’ensemble des mémoires (RAM, ROM, cache, mémoire de masse,...) pour la mémorisation des instructions mettant en œuvre le procédé selon l’invention, des paramètres et réglages utiles à cette mise en œuvre et pour la mémorisation des résultats des calculs intermédiaires et finaux.

Optionnellement, l’appareil photographique 12 est avantageusement logé dans un boîtier de protection 26 muni d’une fenêtre transparente 28 en face de l’objectif 16 pour le protéger des intempéries et ainsi permettre de laisser le dispositif à demeure sur le site choisi. Dans une variante, l’unité 24 est également logée dans le boîtier de protection 26 et communique, via par exemple un réseau sans fil (téléphonique, satellitaire, etc...), le résultat du procédé qu’elle met en œuvre à un site distant où la production du parc photovoltaïque est gérée.

En fonctionnement, un flux d’images numériques horodatées {INn} est donc produit par l’appareil 12, une image numérique INn étant acquise par exemple toutes les Te secondes, l’entier n représentant le nieme instant d’acquisition nxTe depuis un instant initial n = 1, et donc la //eme image dans le flux d’images {INn} . En variante, l’acquisition est réalisée de manière non régulière avec un pas d’acquisition non constant dans le temps.

Chaque image INn comprend A lignes de pixels par B colonnes de pixels et comprend trois canaux d’un espace colorimétrique, par exemple les canaux rouge (r), vert (v) et bleu (b) codés respectivement par trois matrices R, V, B de M'4xB. Les valeurs numériques composant les matrices sont codées sur un nombre bits prédéterminé, par exemple 8 bits, correspondant à 256 niveaux de valeurs (de 0 à 255). Par la suite, les canaux r,v,b d’une image en couleur quelconque « IN » sont respectivement notés /?(/ÏV), V(IN) et B(IN).

Dans ce qui suit l’espace colorimétrique considérée sera l’espace RVB, étant entendu que l’invention s’applique à d’autres espaces colorimétriques, par exemple l’espace H SV (teinte, saturation, intensité lumineuse). De même, il est décrit un nombre de bit de codage égal à 8, correspondant à des valeurs comprises entre 0 et 255. Bien entendu, l’invention s’applique à un nombre de bits quelconque. C. Traitement informatique

Il va à présent être décrit, en relation avec les organigrammes des figures 2A-2E, un procédé selon un mode de réalisation de prévision de la production d’électricité d’un parc photovoltaïque, et plus particulièrement la chute d’ensoleillement sur le parc photovoltaïque entraînant une chute de la production, mis en œuvre par le dispositif de la figure 1, ce procédé se fondant sur l’estimation selon l’invention d’une direction moyenne de déplacement et d’une vitesse moyenne des nuages dans le ciel. C.l. Calibration

Le procédé débute par une phase 40 de positionnement et de calibration du dispositif d’imagerie 10.

Une première étape 400 consiste à charger dans la mémoire de l’unité 24 l’ensemble des instructions et paramètres utilisés par le procédé, notamment : - un modèle astronomique de la position du soleil dans le ciel, par exemple déterminant l’azimut et l’élévation du soleil dans le ciel en fonction de la longitude, la latitude, de l’heure et du jour de l’année. Ce modèle est par exemple le « solar position calculator » de la NOAA (« National Oceanic and Atmospheric Administration ») ; - un modèle de la formation d’image sur le capteur 14 par l’objectif 16. Pour le fisheye, il s’agit par exemple du modèle décrit dans le document « Géométrie modelling and calibration of fisheye lens caméra System », Schwalbe E., In Proc. 2nd Panaramic Photogrammetry Worshop, Int. Archives of Photogrammetry and Remote Sensing (vol. 36, No. Part5, p. W8), février 2005. Ce modèle permet notamment de modéliser l’image d’un objet photographié par le dispositif 12 sur le plan focal de l’objectif 16, et donc sur le capteur 14.

En connaissant la position (θν θ2, Az) du capteur 14 dans un repère terrestre, et donc la position du capteur par rapport au soleil, la date et l’heure, la combinaison du modèle astronomique et du modèle de la formation d’image permet de calculer l’image du centre du soleil sur le plan focal de l’objectif, et donc le centre du disque solaire dans les images produites par l’appareil 12. Cette combinaison est ci-après désignée par les termes « modèle astronomique optique ».

En 402, un opérateur positionne l’appareil 12 sur un site choisi avec l’objectif visant le ciel. La position de l’appareil est ainsi défini par trois paramètres (01; θ2,Αζ), le couple d’angles (θ±, θ2) définissant la position du capteur 14 par rapport à un plan horizontal, et le paramètre Az représentant l’angle que forme l’axe passant par les centres des bords haut et bas du capteur avec un azimut, par exemple l’azimut nord. L’opérateur cherche notamment à positionner le capteur 14 sur un plan horizontal avec l’azimut nord correspondant au milieu du bord supérieur des images produites par l’appareil 12. La position souhaitée, i.e. idéale, de l’appareil est alors égale à (θνθ2,Αζ) = (0,0,0). Les coordonnées géographiques, e.g. la longitude et la latitude, du dispositif sont par ailleurs mémorisées dans l’unité 24.

Evidemment, il est recherché à positionner le dispositif dans une zone dégagée de manière à maximiser la portion de ciel dans les images. Toutefois, le dispositif peut être implanté dans une zone très contrainte en termes d’éléments non pertinents présents dans le champ de l’appareil (ou « objet non météorologiques », ou « onm », par exemple des arbres, des constructions, des pilonnes, etc...), l’invention permettant en effet une prévision robuste y compris sur des portions minimes de ciel et y compris avec le soleil caché.

Une fois le dispositif 10 positionné, une calibration 404 est mise en œuvre. Notamment, le positionnement parfait de l’appareil est difficile à obtenir, de sorte qu’il présente une position réelle (θ[, θ2,Αζτ) différente de la position théorique (θ1,θ2,Αζ). La position du dispositif étant utilisée lors d’étapes de traitement selon l’invention, une correction de l’écart entre les positions réelle et théorique est mise en œuvre. L’étape de calibration 402 a pour objectif de déterminer les paramètres de cette correction et est décrite plus en détail par la suite. La position corrigée (0[, θ2, Azr) est mémorisée par l’unité 24 et est celle utilisée par le modèle astronomique optique pour calculer le centre du disque solaire dans les images produites par l’appareil 12.

La phase initiale 40 se poursuit par le calcul, en 406, d’un masque d’exclusion des objets non météorologiques de la scène présents dans le flux d’images acquis. Ce masque permet de définir les pixels des images sur lesquels portera le traitement pour la prévision de l’ensoleillement. Etant dans une phase d’initialisation et de calibration, il est possible de mettre en œuvre toute méthode connue permettant de déterminer ce type de masque, y compris les méthodes les plus sophistiquées et consommatrices de temps et de ressources matérielles, voire même avec l’aide de l’opérateur qui désigne à la main tout ou partie de ces objets.

De manière avantageuse, la détermination du masque d’exclusion se fonde sur les spécificités des images acquises à savoir celle du ciel, ce qui permet une détermination automatique et rapide du masque. Notamment, en journée, le ciel « clair » (portion de ciel sans nuage et sans le disque solaire, i.e. les portions de ciel bleu), les nuages, le disque solaire et les objets non météorologiques ont des luminosités différentes, comme illustré par exemple à la figure 3.

La détermination du masque d’exclusion comporte ainsi : a) l’acquisition par l’appareil 12 d’une image numérique INm en journée, par exemple par temps ensoleillé ou légèrement nuageux avec le soleil haut dans le ciel, e.g. à son zénith ; b) le calcul d’une matrice de luminosité ILm de l’image INm, d’une manière connue en soi, par exemple selon la formule ILm = - X (Æ + V + B), obtenant ainsi une image en niveau de gris telle qu’illustrée à la figure 3 ; c) la détermination d’un premier masque Mx G RAxB correspondant à l’image optique de la scène sur le capteur, à savoir m^a, b) = 1 si le pixel (a, b) de INm appartient à l’image de la scène sur le capteur 14 (i.e. au disque fisheye sur l’image) et 0 sinon ; d) l’application d’un algorithme de segmentation d’histogramme de luminosité, tel que par exemple décrit dans le document « A threshold sélection method form gray-level histograms», de Otsu N, Automatica, 11(285-296), 23-27, 1975. Ce type d’algorithme consiste à diviser l’histogramme de la matrice Mx 0 ILm en une pluralité de parties homogènes du point de vue de leurs valeurs, par exemple quatre parties en calculant trois seuils SI, S2, S3 tels que: (51,52,53) = argmin (σ(0,$1) + σ($1 + l,s2) + a(s2 + l,s3) + σ($3 + 1,255)) (sl,s2,S3)e[0,255]3 où σ(η1,η2) est l’écart type des valeurs de Γhistogramme d’abscisses comprises entre ni et n2. En adoptant 5 niveaux de gris NG1 0 <NG2 < NG3< NG4<NG5= 255, l’image M1 0 lL°mU résultante de cette segmentation (iZ£ftu(a, b) = N G1 si ilm(a, b) < 51, iZ£ftu(a, b) = NG2 si 51 < ilm < 52, etc.) est illustré à la figure 4. Ce type de problème est simple de résolution à l’aide d’un solveur approprié, par exemple mis en œuvre par l’unité 24 ; e) Pour chacun des angles d’un cercle trigonométrique centré sur le centre de l’image (entre 0° et 360° avec un pas angulaire p prédéterminé, par exemple de 1°), il est déterminé le pixel de la matrice 0 ILm le plus proche du centre du cercle dont la luminosité ZZm(a,b) est inférieure au seuil SI qui définit la classe de luminosité la plus sombre issue de la segmentation. Il est ainsi obtenu un profil, par exemple mémorisé dans un tableau de 360//? valeurs d’angles, de la distance au centre de l’image qui est la plus petite distance identifiée d’un pixel sombre, le profil résultant de l’image Mt 0 ILm étant illustré à la figure 5 ; f) la détermination d’un second masque M2 G M/1xBavec m2(a,b) égal à 0 pour tous les pixels au-delà et sur le profil et égal à 1 sinon, le masque M2 étant illustré à la figure 6.

Le masque M2 est dans une première variante celui utilisé lors du traitement ultérieur de prévision d’ensoleillement. Toutefois à la figure 6, on constate que les limites de la zone masquée (en noir) sont très proches par endroits des cimes des arbres et des contours des branches. On se doute qu’en cas de vent, ce masque ne permettrait pas d’exclure les branches qui se déplaceraient et apparaîtraient alors dans la partie visible (en blanc) du masque. Dans une seconde variante, une correction supplémentaire du masque M2 est appliquée pour tenir compte de ce phénomène. Notamment, un algorithme de dilatation est mis en œuvre sur le masque M2 avec un rayon suffisant pour englober au moins le feuillage, par exemple au moyen d’un masque de dilatation circulaire ayant un rayon compris entre 10 et 50 pixels. Plus particulièrement, pour chacun des pixels préalablement défini comme masqué du masque M2 (m2(a, b) = 0), un masque circulaire est centrée sur ce pixel et tous autres pixels de M2 contenu dans le masque circulaire sont réglés sur 0. Le masque M3 résultant de la dilation du masque M2 à l’aide d’un masque de dilation circulaire de rayon égal par exemple à 10 pixels est illustré à la figure 7. En se référant à la représentation de la matrice M3 0 ILm de la figure 8, on peut ainsi constater que le masque d’exclusion permet d’identifier les parties utiles de l’image (ici le ciel, le nuage et le soleil) sur lesquels porteront les traitements ultérieurs mis en œuvre dans le cadre de la prévision de la chute d’ensoleillement sur le parc photovoltaïque.

Le masque d’exclusion correspondant aux objets non météorologiques Monrn = M2 ou Monm = M3 est alors stockée dans une mémoire de l’unité 24.

Une fois la phase initiale de positionnement et calibration terminée, le dispositif est mis en service. L’appareil 12 produit le flux d’images numériques {INn} qui est traité en 50 par l’unité d’information 24 pour l’estimation de l’ensoleillement. Le procédé ayant pour but de déterminer les variations infra-journalières d’ensoleillement, le traitement complet est mis en œuvre avantageusement uniquement lorsqu’il existe du ciel clair. Par exemple, l’unité 24 détermine la luminosité moyenne de chaque image INn et l’unité 24 déclenche le traitement si la luminosité moyenne (e.g. de l’image Monm 0 ILm) dépasse un seuil prédéterminé, signifiant que le ciel n’est pas couvert. D’autres tests, ou aucun test, de déclenchement peuvent bien entendu être mis en œuvre. L’étape de traitement informatique 50 comprend, pour chaque nouvelle image numérique acquise INn, une étape 500 de prétraitement de l’image, suivie d’une étape 502 d’estimation de la direction moyenne des nuages à l’instant d’acquisition n de l’image, suivi d’une étape 504 d’estimation de la vitesse moyenne des nuages audit instant n, suivi par une étape 506 de prévision de la chute de l’ensoleillement au dessus du système 12, et donc du parc où ce dernier est implanté, en fonction de la direction et de la vitesse moyenne calculée à l’instant d’acquisition n.

C.2. Prétraitement de l’image INV

Afin de faciliter la segmentation ultérieure entre ciel clair et nuages, le procédé débute, en 500-2, par l’exclusion supplémentaire du disque solaire si ce dernier est présent au moins partiellement dans l’image INn nouvelle acquise, ce qui ne laisse donc subsister que les portions de l’image acquise INn correspondant au ciel clair et aux nuages. Plus particulièrement, l’unité 24 calcule la position théorique Pn du centre du disque solaire dans le plan focal et donc le cas échéant dans l’image INn en fonction du modèle astronomique optique, puis détermine un disque de rayon rs ayant pour centre la position Pn. Le rayon rs, paramètre enregistré dans l’unité 24, est une valeur déterminée par exemple lors d’une campagne de test antérieure et est supérieure ou égale au plus grand rayon de disque solaire observé au cours de l’année, ou d’une période définie de l’année, sur toute l’année ou sur plusieurs années. D’autres méthodes de calcul dudit rayon sont évidemment possibles. Le rayon rs permet donc à d’exclure de manière certaine l’image du soleil dans l’image numérique INn sans pour autant exclure une zone trop importante. Cette valeur, qui dépend de la résolution de l’image, est par exemple égale à 80 pixels.

Un masque d’exclusion M4n de RAxB est alors calculé par l’unité 24. Ce masque est égal au masque Monm lorsque le disque solaire n’est pas compris dans l’image. Dans le cas contraire, le nouveau masque d’exclusion M4n correspond à l’exclusion des objets non météorologiques et des pixels du disque solaire compris dans l’image INn, à savoir m4n(a, b) = 1 si m3(a, b) = 1 ou si (a, b) appartient au disque solaire, et m4n(a, b) = 0 sinon. L’image M4n 0 ILn correspond ainsi uniquement aux portions de ciel clair et de nuages de la scène photographiée et est illustrée à la figure 9.

Une segmentation de l’image M4n 0 INn est ensuite mise en œuvre, en 500-4, pour identifier dans celle-ci les portions de ciel clair et les nuages.

Avantageusement, la matrice M4 n 0 INB~B est utilisée pour cette segmentation, où IN*~B = |R(/Nn ) — β(/Αη )|. Comme décrit ci-dessous en relation avec la calibration de la position réelle de l’appareil 12, la matrice IN£~B, issue de la différence des canaux rouge et bleu, fait apparaître de manière contrastée les portions de ciel clair et les nuages. Un simple seuillage sur la valeur centrale de l’histogramme de cette matrice peut ainsi suffire à segmenter l’image. Un algorithme de type Otsu tel que décrit précédemment, avec un unique seuil, est par exemple alternativement mis en œuvre pour distinguer dans l’histogramme la classe du ciel clair et la classe des nuages. Ceci permet une segmentation plus précise, notamment concernant les portions de nuages. Suite à cette segmentation, une matrice binaire M4 n 0 IN£cn est ainsi calculée (« CCN » pour « ciel clair/nuage ») avec par exemple in^CN (a, b) = 1 si le pixel (a, b) appartient à un nuage et 0 si le pixel (a, b) appartient à une portion de ciel clair. Cette matrice est illustrée à la figure 10.

Le prétraitement 500 se poursuit en 500-6 par la correction de la distorsion géométrique induite par l’objectif 16, notamment la correction du masques M4n et de l’image binaire INfiCL. Ce type de traitement, également connu sous le terme de « redressement », vise à ce qu’une droite dans la scène photographiée soit également une droite dans l’image redressée.

De manière privilégiée, le redressement mis en œuvre se fonde sur le fait que l’invention est peu sensible aux effets de perspectives dans l’image, et utilise donc avantageusement le redressement simplifié de la projection d’une scène 3D sur un plan 2D par un fisheye, par exemple celui décrit dans le document de Schwalbe E. cité précédemment Ce redressement, qui ne corrige que les déformations de l’image liées à la projection au travers de l’objectif, est donc très simple et rapide à calculer.

Dans ce qui suit, les portions de ciel clair n’étant plus utilisées en tant que telles pour la détermination de la direction des nuages, le procédé se poursuit par exemple par la génération d’une image binaire redressée, notée IN£, identifiant les nuages de l’image INn, à savoir in£(a, b) = 1 si le pixel (a, b) appartient à un nuage et 0 sinon. Cette image inclut donc le masque M4 n.

Enfin, l’unité 24 calcule en 500-8 la position Sn = (xn,yn) du centre du disque solaire dans l’image redressée IN£. Le prétraitement mis en œuvre par l’unité 24 calcule ainsi pour chaque image numérique INn acquise par le dispositif 12, une image binaire redressée IN%, ainsi que la position Sn du centre du disque solaire dans l’image redressée.

On notera que les dimensions des images redressées ne sont pas nécessairement égales à celles des images issues du capteur. Pour des raisons de concision, il sera par la suite fait cependant cette hypothèse, les matrices et calculs pouvant bien évidemment être définis pour des dimensions différentes. La même hypothèse est utilisée pour la rotation d’une matrice.

Le procédé continue par l’estimation de la course moyenne des nuages à l’instant n, et plus particulièrement par l’estimation 502 de leur direction moyenne, suivi de l’estimation 504 de leur vitesse moyenne selon la direction moyenne. C.3. Estimation de la direction moyenne des nuages L’estimation de la direction moyenne des nuages consiste dans une première étape à découper la nouvelle image IN£ issue du prétraitement de l’image INn en une pluralité de zones, chacune associée à une direction particulière de recherche Dk prédéterminée et constante dans le flux d’images acquises {INn} . Puis, dans une seconde étape, l’estimation détermine l’évolution des nuages en fonction du temps dans chacune des zones. La direction Dk de la zone de variation nuageuse minimale correspond alors à la direction moyenne des nuages. De manière avantageuse, l’estimation selon l’invention est configurée pour tenir compte de la prévision d’ensoleillement recherchée.

En se référant à la figure 11 illustrant un rectangle sur une image M4n®/Nn redressée, l’étape 502 d’estimation de la direction moyenne des nuages débute, en 502-2, par la définition d’un ensemble de rectangles de détection {/?fcn} de longueur et de largeur en pixels identiques mais d’orientations différentes. Plus particulièrement, les rectangles de détection Rkn ont : - une largeur en pixels / égale au moins égale au diamètre du halo solaire, noté D%. Le halo solaire, également appelé « cône solaire », correspond à une zone de luminosité très proche de celle du disque solaire et entourant ce dernier, dont l’occultation par des nuages induit donc également une chute de l’ensoleillement. La largeur / est par exemple déterminée lors d’une campagne de test antérieure, mémorisée dans l’unité 24 et est supérieure ou égale au diamètre maximal du halo solaire observé lors de cette campagne. Par exemple, pour une image IN£ de 500 pixels par 500 pixels, la largeur / est égale à 160 pixels ; - un petit côté centré sur le centre Sn du disque solaire. Le rectangle permet ainsi de tenir compte de l’ensemble des nuages susceptibles d’intercepter ultérieurement le halo solaire D%, et ainsi faire varier l’ensoleillement ; - une longueur en pixels L, par exemple égale à la plus grande dimension des images IN£ (e.g. 500 pixels dans les exemples illustrés sur les figures). Comme détaillé ci-après, la longueur des rectangles définit la fenêtre temporelle de prévision des variations de l’ensoleillement. Notamment, plus la longueur L est importante et plus les nuages éloignés du soleil sont pris en compte pour la prévision d’ensoleillement, notamment la chute d’ensoleillement; - une direction prédéfinie Dk, parallèle à la longueur L, repéré par un angle ak par rapport à un axe prédéterminé de l’image, par exemple un axe horizontal des images (e.g. leur bord inférieur). Les directions Dk sont par exemple régulièrement espacées angulairement entre 0° et 360°, avec un pas angulaire Θ qui règle la précision souhaitée sur la direction moyenne recherchée. Par exemple, le pas Θ est choisi égal à 1°, résultant en K = 360 directions Dk fixes étudiées.

Comme illustré aux figures 12A-C, en fonction de la position du disque solaire dans l’image redressée, un rectangle peut être compris entièrement dans l’image (figure 12A), ou partiellement dans l’image si le disque solaire est dans l’image et proche d’un bord de celle-ci (figure 12B) ou si le disque solaire est hors de l’image et proche d’un bord de l’image (figure 12C).

Le procédé consiste, dans une étape suivante 502-4, à déterminer, pour chaque zone Zk)îl de l’image binaire redressée IN£, comprise dans un rectangle Rkn, la distribution DDkn des pixels (DD pour « distribution déplacement ») correspondant aux nuages le long de la largeur / du rectangle Rk n, et donc dans une direction normale à la direction Dk associée à ce rectangle.

Par exemple, pour chaque rectangle de détection Rk n, le procédé : produit une image binaire rectangulaire lRkiU correspondant au contenu du rectangle Rkn (i.e. la zone Zfcn, le cas échéant complétée d’une portion de pixels à 0 correspondant à la zone du rectangle Rkn non comprise dans l’image IN£, tels qu’illustré aux figures 12B et 12C) ; applique une rotation d’angle 90° — ak afin d’obtenir une image rectangulaire IRkn de l x L pixels, le petit côté du rectangle étant ainsi horizontal ; calcule la distribution DDkn. Par exemple, la distribution DDkn est un histogramme comptabilisant pour chaque colonne de l’image lRkn le nombre de pixels de la colonne qui sont égaux à 1, i.e. correspondant à des nuages.

Pour chaque image numérique INn acquise par le dispositif 10, l’unité 24 calcule donc un ensemble de distributions {DDkn} respectivement associées à des directions fixes {Dk} . L’unité 24 produit donc, pour le flux d’images {INn} , K flux de distributions, par exemple noté pour le £ieme flux FDDkn = {DDkl, DDk2, —, DDkn_1,DDkn}. Des exemples de flux de distribution FDDkn sont illustrés graphiquement aux figures 13A-13J qui correspondent respectivement aux angles ak égaux à 0°, 10°, 20°, 30°, 40°, 50°, 60°, 70°, 80° et 90°. En se référant à la figure 14, la représentation graphique consiste à transformer l’échelle des valeurs des distributions (de 0 à I concernant des histogrammes) en échelle de gris, ce qui donne une représentation graphique d’une distribution sous la forme d’une étiquette, puis à juxtaposer les étiquettes en fonction des valeurs croissantes de n.

Selon l’invention, l’analyse de la variabilité des flux de distributions FDDkn permet de déterminer la direction moyenne des nuages à l’instant n. A cet effet, dans une étape suivante 502-6 du procédé, la variabilité à l’instant n est calculée en fonction des p dernières distributions calculées {DDkn_p, DDfcn_p+1, ...,DDkn_v DDkn], i.e. sur une fenêtre temporelle [(n — p) xTe;nx Te\, où p est entier prédéterminé. La durée de la fenêtre est par exemple comprise entre 4 secondes et 44 secondes, ce qui correspond à p= 1, 2,.., ou 11 pour un pas d’acquisition Te de 4 secondes. La variabilité est par exemple calculée par tout critère approprié CVkn, notamment un critère de ressemblance ou de dissemblance entre les distributions comprises entre l’instant d’acquisition n et l’instant antérieur d’acquisition n-p.

Avantageusement, le critère CVkn est égal à la somme des ressemblances entre deux distributions espacées de m > 1 pas de temps, à savoir un critère de dissemblance CVkn pour m = 1 selon la relation :

(1)

Les inventeurs ont découvert qu’un tel critère comprend deux optimums (deux maxima si critère de ressemblance ou deux minima si critère de dissemblance) en fonction de l’indice k. Ces angles correspondent respectivement à la direction l’angle des nuages approchant le disque solaire et l’angle des pixels nuages s’éloignant du disque solaire (i.e. correspondant à un angle ak et un angle ak + 180°). Autrement dit si les nuages viennent du nord-ouest, la ressemblance des distributions est maximale pour un angle de détection akn de 135° mais également pour l’angle 135°+180°= 215° (direction sud-est) puisque, selon cette direction, les mêmes nuages qui, une fois passés sous le disque solaire, s’en éloignent avec la même direction.

Le procédé détermine ainsi l’indice k, correspondant à l’angle ak, qui maximise (resp. minimise) le critère CVkn entre 0° et 180°, ainsi que l’indice k\ correspondant à l’angle ak + 180°, qui maximise (resp. minimise) le critère CVkn entre 180° et 360°. Par exemple, pour un pas d’angle Θ de 1°, le procédé détermine l’indice k=kopt compris entre 0 et 180 qui maximise (resp. minimise) le critère CVkn = CVkn + CVk+18Q n. La direction Dkoptn correspondant à l’indice kopt déterminé par l’unité 24 pour l’instant n, est ainsi la direction moyenne de déplacement des nuages à l’instant n.

La figure 15 illustre des critères CVkn se fondant sur la relation (1) pour différentes écartements temporels entre les distributions. Les ronds indiquent la valeur de l’angle (Xkopt Pour laquelle le critère CVkn est minimal, et donc la direction Dkoptn de propagation des nuages à l’instant n finalement identifié par le procédé selon l’invention. On constate que le minimum de ces sommes est faiblement sensible à l’écart temporel entre deux distributions. Notamment, l’angle de propagation des nuages varie entre 30 et 28 degrés quel que soit la durée de la fenêtre. D’autres critères CVkn que celui de la relation (1) peuvent bien entendu être utilisés. Notamment d’autres types de norme que la norme L1 d’un vecteur d’un vecteur peuvent être utilisées. De même, en variante un critère R2 ou de Nash est utilisé. Les résultats de ces critères de ressemblances sont illustrés respectivement aux figures 16 et 17. Tous deux sont d’autant meilleurs qu’ils approchent la valeur de 1, le procédé recherchant cette fois un maximum et non un minimum comme précédemment. On constate que ces deux critères ont un maximum peu marqué, ce qui peut rendre sa détermination plus difficile. Le critère selon la relation (1), qui présente un minimum très marqué, est donc privilégié.

La figure 18 illustre une direction de propagation déterminée manuellement par un opérateur. Sur deux images espacées de 7 pas de temps (soit 28 secondes) l’opérateur a dessiné à la main trois vecteurs de déplacements de nuages. L’opérateur a identifié dans les deux images deux zones correspondant aux mêmes nuages. L’origine d’un vecteur correspond ainsi à un pixel de nuage dans la première image et l’extrémité du vecteur correspond au pixel correspondant dans la seconde image. Ces trois mesures de direction de déplacement des nuages sont illustrées par les trois rapporteurs et donnent des valeurs sensiblement identiques et cohérentes avec la direction déterminée par le procédé selon l’invention (30°), étant entendu que le positionnement des vecteurs par l’opérateur est incertain à quelques pixels près en raison de la déformation des nuages.

La direction moyenne de propagation des nuages Dkoptn étant ainsi connue à chaque instant n, il est possible de mettre en œuvre une détermination de leur vitesse moyenne, selon tout traitement connu de l’état de la technique. Un calcul selon l’invention de cette vitesse est à présent décrit. C.4. Estimation de la vitesse moyenne de déplacement des nuages

Une fois déterminée la direction moyenne de déplacement Dfloptn des nuages à l’instant n, avantageusement, mais optionnellement, selon le procédé décrit ci-avant, le procédé se poursuit par l’estimation 504 de leur vitesse moyenne selon ladite direction et à cet instant.

Plus particulièrement, cette estimation débute par une étape 504-2 consistant, pour image binaire redressée IN£, issue du prétraitement de l’image nouvellement acquise INn, à déterminer la distribution DVkoptn (DV pour « distribution de vitesse ») des pixels de la zone Zkoptn,n (i-e· la zone de l’image IN£ associée à la direction Dkoptn) correspondant aux nuages, le long de la longueur L du rectangle /?fcop(. n, et donc selon la direction DkoptiU. Ce calcul est par exemple réalisé sur l’image binaire rectangulaire IRkoptin après rotation de l’angle 90° — akopt en calculant un histogramme comptabilisant, pour chaque ligne de l’image IRkopt/n après rotation, le nombre de pixels de la ligne qui sont égaux à 1, et donc un histogramme selon une direction parallèle à la direction Dk n, tel qu’illustré à la figure 19.

Un flux de distributions FDVkopt§n = [DVkoptil, DVkopti2, -,DVkoptin_ltDVkopt;n^ est ainsi obtenu par l’unité 24, tel qu’illustré à la figure 20 en adoptant la même représentation graphique que précédemment avec l’axe des abscisses correspondant aux instants n d’acquisition. Comme illustré sur cette figure, on constate que ce flux comprend des structures affines qui apparaissent par intermittence. Ces structures correspondent à des passages de nuages. Elles apparaissent d’abord pour des ordonnées hautes (i.e. un banc de nuage apparaît dans le haut de la vignette) puis décroissent en ordonnée au fur et à mesure que ce banc de nuage se rapproche du disque solaire (le centre du disque solaire ayant la coordonnée 0), c’est à dire diminue en ordonnée au fil des images successives dans le flux d’images {INP}. Les inventeurs ont constaté que la pente de ces structures est directement liée à la vitesse moyenne V de déplacement des nuages selon la direction Dkop( n (en pixels par seconde) selon la relation :

(2) où β est l’angle que forme la direction des structures affines avec l’axe des abscisses.

Dans une étape 504-4 suivante, le procédé détermine donc la valeur de l’angle β des structures affines en fonction des q dernières distributions calculées {DVkopt:n-p,DVkoptin_p+1,...,DVkoptin_1,DVkoptin], i.e. sur une fenêtre temporelle [(n — q) x Te; η X Te], où q est entier prédéterminé. La durée de cette fenêtre temporelle est déterminée, par exemple lors d’une campagne de test précédent. Par exemple, une fenêtre temporelle comprise entre 4 secondes et 44 secondes est choisie, ce qui correspond à q= 1, 2,..,11 pour un pas d’acquisition Te de 4 secondes, ce qui a priori est suffisant pour une détermination avec suffisamment de précision de la vitesse moyenne, tout en s’assurant que la vitesse des nuages ne subissent pas de grande modification. D’autres durées sont possibles, par exemple plus importantes. Tout algorithme approprié mis en œuvre par l’unité 24 peut convenir à cet effet (par exemple une régression linéaire).

On note cependant sur la figure 20 que la base des structures (zones proches de la coordonnée 0) peut être plus ou moins large en raison de la proximité avec le soleil. Cet évasement peut perturber la détermination de l’angle par des algorithmes classiques. De manière privilégiée, une approche similaire à celle utilisée pour l’identification de la direction Dkoptn est mise en œuvre. Plus particulièrement, différentes rotations du flux de distribution FDVk n sont testées (correspondant graphiquement à differentes rotations des axes de la figure 20). Normalement, pour une certaine rotation de ces axes, les structures affines deviennent verticales. Dans ce cas, la distribution sur l’axe des abscisses des pixels des nuages présente des pics de valeurs supérieures aux pics des distributions sur l’axe des abscisses des pixels nuages de rotations différentes (les pics sont les plus faibles quand les axes seront tournés d’une valeur conduisant à des structures affines parallèles à l’axe des abscisse). A cet effet, le procédé selon l’invention comprend par exemple : - la production d’une matrice DVn = [DVkoptn_q DVkoptq+1... DV^^DV^^], les distribution DVkoptJl étant codées sous la forme de vecteur colonne de dimension L l’application de différentes rotations d’angles βη prédéterminés et mémorisés dans l’unité 24, et de centre par exemple de coordonnées (n-q, L) à la matrice DVn. Les angles βη sont espacés angulairement entre 0° et 90°, avec un pas angulaire qui règle la précision souhaitée sur la vitesse moyenne recherchée, par exemple un pas égal à 0,5°; le calcul d’un histogramme Histnu pour chaque matrice DVnu obtenue par rotation et le calcul du maximum Mnu de cet histogramme, à savoir le maximum du vecteur [Σ)=ι dVuO-’j') Σ)=ι dvu(2,j) ... Σ)=ι dvu(n-q,j)]; la détermination de l’angle βηορ(,η comme étant l’angle βη correspondant au maximum maxu Mn udes Mn u calculés.

Le procédé calcule alors en 504-6 la vitesse moyenne Vn des nuages, en pixels par seconde, en se fondant sur la relation (2).

La portion supérieure de la figure 21 illustre graphiquement la rotation de la matrice DVn illustrée à la figure 20 pour un angle βη. La figure 22 illustre l’évolution du maximum Mnu en fonction de l’angle βη. On constate que le maximum de la courbe Mnu(/?u) présente un maximum marqué correspondant à l’angle βηορί,η· La figure 23 illustre la détermination manuelle par un opérateur de l’angle βηορΙ,η qui trace les pentes des principales structures de la figure 20. Le rapporteur à droite de la figure 23 montre que l’angle /?optnobtenu par le procédé selon l’invention est cohérent avec une mesure de celui-ci faite directement sur l’image des distributions DVk n.

En outre, l’écart de l’angle /?Uop(.n par rapport à un angle de référence (e.g. 180° si les angles sont mesurés entre 0° et 360° ou 0° si les angles sont mesurés entre -180° et 180°) donne le sens de propagation des nuages selon la direction Dk n. C.5. Prévision de l’ensoleillement

Une fois déterminée la vitesse moyenne des nuages à l’instant n, le procédé se poursuit, en 506, par l’estimation de l’ensoleillement, en particulier l’estimation de la chute d’ensoleillement en raison de nuages passant dans le cône solaire.

En se référant à la figure 24 qui illustre une partie du flux FDVk n, l’axe des ordonnées représente la distance (en pixels) à la base des rectangles de détection (^/copt,n} , c’est-à- dire la distance (en pixels) au centre du disque solaire dans le flux d’images {INn} . Cette distance peut être transformée en temps puisque la vitesse moyenne de déplacement Vn des nuages est connue. Par exemple, en reprenant la vitesse calculée en fonction de la figure 20 (1,48 pixel/seconde), un déplacement de 100 pixels sur l’axe des ordonnées correspond à une durée d’environ 100/1,48 ~ 68 secondes (soit environ 17 images dans le cas particulier d’une période d’acquisition Te de 4 secondes). Autrement dit, une détection de nuage effectuée à l’ordonnée égale à 100 laisse penser que le nuage interceptera le disque solaire 68 secondes plus tard (ou 17 images plus tard). Cela signifie que dans le cas de cette séquence, l’anticipation temporelle maximale est de 500 / 1,48 ~ 340 secondes, soit 5 minutes 40 secondes (500 étant la taille des images redressées et la longueur Z des rectangles de détection).

Toujours en se référant à la figure 24, on peut constater que les nuages en cours de progression vers le disque solaire correspondent aux valeurs significativement non nulles de la distribution DVk n de la dernière image acquise INn. Connaissant donc la position en pixels des nuages à cet instant (i.e. l’ordonnée des valeurs non nulles de la distribution DVkoptlnl ü est ainsi possible de prédire leur moment d’arrivée dans le cône solaire et sur le disque solaire. Le procédé selon l’invention estime donc le passage de nuages devant le soleil en fonction des ordonnées des valeurs non nulles de la dernière distribution DVkopt n et de la vitesse moyenne Vn de déplacement des nuages.

Par exemple, l’estimation de l’ensoleillement consiste à détecter, en 506-2, les valeurs significativement non nulles de la dernière distribution DVk n et à mémoriser leurs ordonnées respectives {Ys} , puis, à chacune de ces ordonnées Ys, à déterminer, en 506-4,

Y l’instant d’arrivée Ts = — des nuages correspondant dans le halo solaire. Les instants

Vn calculés {Ts} correspondent ainsi aux instants pour lesquels il est prédit une chute de 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 inconnues est détectée. Connaissant les instants d’arrivées des nuages, cette alarme peut donc être désactivée ou bien classée avec un niveau d’alarme inférieur (e.g. non critique).

Pour illustrer la capacité prédictive du procédé, il a été extrait de l’image de la figure 20 les chroniques du nombre de pixels nuages par ordonnée fixe (c’est à dire correspondant à une échéance temporelle d’arrivée à la position du soleil fixe), à savoir 50, 100, 150, 250, 300 et 350 pixels, tel qu’illustré respectivement aux figures 25A-G. La vitesse Vn étant de 1,48 pixel par seconde, soit 5,67 pixels par pas de temps Te, en décalant temporellement les figures 25A-G de manière correspondante (de 50 / 5,67 = 8,81 pas de temps pour la Figure 25 A, de 100 / 5,67 = 17,6 pas de temps pour la Figure 25B, de 150 / 5,67 = 26,5 pas de temps pour la Figure 25C, de 200 / 5,67 = 35,3 pas de temps pour la Figure 25D, de 250 / 5,67 = 44,0 pas de temps pour la Figure 25E, de 300 / 5,67 = 52,9 pas de temps pour la figure 25F, et 350 / 5,67 = 61,7 pas de temps pour la Figure 25G), les pics des figures ainsi décalées doivent se correspondre et correspondre également à des pics de chute d’ensoleillement. La figure 26 correspond à la superposition des figures décalées. Comme on peut le constater, les pics des différentes figures décalées se superposent, vérifiant en outre que la vitesse des nuages est restée sensiblement constante. La figure 27 illustre un signal enregistré par un pyromètre disposé à proximité du système à des fins de validation du procédé selon l’invention lors de l’enregistrement des images ayant données lieu à la figure 20 ainsi que la superposition des figures décalées. Comme cela est connu en soi, un pyranomètre permet de mesurer la quantité d’énergie solaire incidente sur celui-ci. Comme on peut le constater, les chutes de signal du pyranomètre, qui correspondent à des chutes d’ensoleillement, coïncident avec les pics des figures décalées, prouvant ainsi la capacité du procédé à prévoir l’ensoleillement, et en particulier les chutes d’ensoleillement.

D. CALIBRATION DE LA POSITION REELLE DU DISPOSITIF D’IMAGERIE

Comme décrit précédemment, le positionnement optimal de l’appareil photographique est difficile à obtenir, de sorte qu’il présente une position réelle (0[, 0j,i4zr) différente de la position théorique souhaitée (θ1,θ2,Αζ). Il va à présent être décrit l’étape de calibration 402 mis en œuvre par l’unité 24, ayant pour objectif de déterminer les paramètres de cette correction. Cette calibration peut être réalisée au jour le jour, de manière automatique et/ou sur commande d’un utilisateur, après mouvement du dispositif de prise de vue par exemple (suite à un coup de vent notamment). Cette calibration permet également de détecter si l’appareil a bougé. Notamment, en mettant en œuvre la calibration et en comparant la position réelle obtenue avec une position réelle obtenue au cours d’une calibration antérieure, il est déterminé que l’appareil a bougé si la différence entre les deux positions est supérieure à un seuil. En cas de différence importante, une alarme peut être déclenchée pour obtenir le déplacement d’un opérateur afin de remettre d’aplomb le dispositif de prise de vue.

Dans une première étape 402-2, un flux d’images numériques de calibration {INfi} 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 Pfi du centre du disque solaire dans chaque image INfi est calculée en fonction de l’image elle-même. Notamment, le calcul par l’unité 24 de la position Pfi du centre du disque solaire met en œuvre un traitement rapide se fondant sur la spécificité de l’image acquise, à savoir celle du ciel, et donc se fondant sur la présence de zones à très forte dominante bleue (le ciel clair) et de zones ayant des tons plus gris (nuage et soleil). Pour chaque image du flux {INfi} , le calcul comprend tout d’abord le calcul de la matrice de luminosité ILcn et de la matrice INfi'R~B = \R(INfi) — B (INfi) |. Cette dernière permet de séparer au premier ordre le ciel clair du reste de l’image, les fortes différences correspondant en effet aux parties de ciel clair alors que les faibles différences appartiennent aux nuages et au disque solaire. Un balayage pixel par pixel des matrices Monrn 0 INfi’R~B et Monm 0 ILcn est ensuite effectué et : si la luminosité du pixel est suffisamment importante (i.e. iZ£(a, b) > su où su est un seuil proche du maximum max{ Monm 0 ILcn) de la matrice Monm 0 ILcn, avec sa compris par exemple entre 80% et 90% de max( Monm 0 ILcn) ou par exemple égal à 250, et que la valeur absolue de r-b du pixel est suffisamment petite, signifiant que le pixel est plus gris que bleu (i.e. infiR~B(a, b) < sr_b où sr_b est un seuil proche de zéro, avec sr_b inférieur ou égal par exemple à 5), le pixel (a, b) est déclaré comme candidat au disque solaire. Les pixels candidats constituent ainsi un sous-ensemble de l’image Monm 0 ILfi-, - Ce sous-ensemble est érodé pour éliminer les groupements de pixels candidats de trop petites tailles et pour favoriser l’identification ultérieure de zones distinctes de pixels candidats. Un filtre circulaire d’érosion de diamètre de sept pixels est par exemple mis en œuvre. Il est ainsi obtenu un ensemble de zones de pixels candidats ; le centre de gravité et la distance du centre de gravité à une position prédéterminée a priori (par exemple la position du centre du disque solaire calculée à l’image précédente ou la position théorique Pf du centre du disque solaire calculé par le modèle astronomique optique notamment) sont ensuite calculés pour chaque zone de pixels candidats ; - la zone dont le centre de gravité est le plus proche de la position prédéterminée est retenue comme correspondant au disque solaire ; et le plus grand cercle inscrit dans la zone retenue est alors calculé, et correspond au disque solaire dans l’image. La position finale réelle du centre du soleil P£ retenue pour les calculs ultérieurs mis en œuvre par le procédé est alors le centre de ce cercle.

La figure 28 illustre un exemple de positions réelles P£ déterminée pour le flux de calibration {INf} .

Le procédé se poursuit, en 402-6, pour chaque image du flux {INf} , par le calcul de la position théorique Pf(01,02,Az) du centre du disque solaire calculée par le modèle astronomique optique, cette dernière est calculée à l’aide du modèle astronomique optique connaissant l’horodatage des images et en fixant arbitrairement les valeurs des trois angles (θνθ2,Αζ). En testant un nombre suffisant de combinaisons de valeurs de ces trois angles (θ±, θ2,Αζ), il est alors possible de déterminer la combinaison qui produit la meilleure comparaison entre les coordonnées Pf calculées du centre du disque solaire par cette méthode analytique et les coordonnées Pf du centre du disque solaire par calculées lors de l’étape 402-4. Par exemple, la valeur (θ[,02,Azr) des angles finalement retenue est celle minimisant la somme des différences absolues entre l’ensemble des coordonnées {Pf} et {Pf}. D’une manière générale, toute méthode d’optimisation, e.g. mise en œuvre par ordinateur, pour résoudre un problème d’optimisation selon la relation suivant convient :

où Ω est un ensemble de recherche, par exemple un hypercube centré sur la position souhaitée par l’opérateur lors de l’installation du dispositif d’imagerie (e.g. (0,0,0)) ou R3, Vi = 1, ...,n, Pf est un estimateur de la position du disque solaire dans l’image INf se fondant sur une analyse de l’image INf (e.g. la position Pf décrite ci-dessus), et || || est une fonction de coût quantifiant un écart entre deux vecteurs de R3 (e.g. la norme Ll, L2, infinie, etc.).

La figure 29 illustre le résultat de cette calibration en indiquant par une flèche dans l’image l’orientation de l’azimut Nord (également représenté sur la boussole en bas à gauche), ainsi qu’un niveau à bulle symbolisant par un point l’écart à l’horizontalité du dispositif. Sur cette figure, la courbe à l’arrière-plan correspond aux positions du soleil Pn et la courbe au premier plan correspond aux positions du soleil P£ avec les trois angles égaux à (0£, >Azr). On constate que les deux angles définissant les écarts à l’horizontal du fish-eye sont inférieurs à 10°, c’est à dire des angles en effet très petits, compte tenu de la précision découlant de l’expérience (appareil photographique monté sur un pied prenant appui sur revêtement en gravier). Comme on l’observe, les positions théoriques et calculées précédemment empiriquement coïncident très bien. Cette façon de procéder peut être jugée bien plus précise que des tentatives de positionnement à la main, d’un appareil photographique par exemple, sur le terrain et finalement plus simple et plus précise que la méthode classique s’appuyant usuellement sur des images numériques de damiers réguliers pour caler un modèle de distorsion empirique couplé à une projection de perspective.

Il a été décrit un mode de réalisation privilégié dans lequel la direction et la vitesse moyenne des nuages sont calculées pour chaque image acquise, ou de manière équivalente à chaque instant d’acquisition. En variante, un flux d’images numériques est acquis sur une durée prédéterminée jugée appropriée pour le calcul (e.g. d’une durée comprise entre 4 secondes et 44 secondes avec un nombre minimal d’images, compris par exemple entre 2 et 11) et mémorisé dans l’unité de traitement. Une fois la dernière image du flux acquise, l’unité calcule la direction et la vitesse comme décrit précédemment (découpage en zones, calcul de distributions, analyse de la variabilité sur le flux acquis, etc...) en fonction de toutes les distributions calculées sur l’ensemble de ce flux. Le procédé selon l’invention est par exemple mis en œuvre régulièrement (par exemple toutes les 5 minutes ou à chaque fenêtre temporelle acquise) ou sur demande d’un utilisateur. Ceci permet par exemple de libérer les ressources de l’unité de traitement pour d’autres types de tâches.

De même, des traitements supplémentaires peuvent être appliqués. Par exemple, une détection de valeurs aberrantes peut être mise en œuvre (valeurs induites par exemple par le passage d’un oiseau, d’une feuille, etc...) et les images ou les distributions correspondantes écartées. Le calcul est alors par exemple effectué sur une fenêtre temporelle plus importante pour conserver le même nombre d’images ou sur la même fenêtre temporelle mais avec un nombre de distributions moins important.

Claims (16)

1. REVENDICATIONS

   1. Procédé d’estimation de la vitesse moyenne Vn de déplacement d’objets en mouvement selon une direction Dkoptn 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 Dkoptn en fonction du flux d’images numériques {INn} ; pour chaque image lNn du flux d’images numériques {INn} : o le calcul d’une image segmentée IN% identifiant les objets en mouvement dans l’image lNn; o le calcul d’une distribution DVkoptJl des valeurs des pixels d’au moins Une zone Zk n de l’image segmentée IN%, la distribution DVkopt>n étant calculée le long de la direction moyenne Dkgptn déterminée ; le calcul de l’inclinaison d’au moins une structure affine comprise dans une matrice DVn = [DVkgptJl_qDVk opt.n-q+1 — DVfcopt.«-lDVfcopt*n] ÎSSUe de 1&amp; concaténation des distributions DVkgptn des images entre la (n-q)ieme image du flux jusqu’à la nlè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.
2. 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 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 {/fVn} .
3. 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 βη prédéterminés à la matrice DVn; le calcul d’un histogramme Histnu pour chaque matrice DVnu obtenue par rotation et le calcul du maximum Mn u de cet histogramme ; la détermination d’inclinaison comme étant égal à l’angle Pu0pt,n correspondant au maximum des maxima Mn u calculés.
4. 4. Procédé selon l’une quelconque des revendications précédentes, caractérisé en ce que l’estimation de la direction Dk 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 {INn} : o le calcul d’un ensemble de zones d’image {Zkn} respectivement associées à des directions d’analyse {Dk} formant des angles {afc} différents prédéterminés avec un axe de référence Aref des images numériques: o pour chaque zone Zkn dudit ensemble de zones {Zfen} , le calcul d’une distribution DDk n des valeurs des pixels de ladite image segmentée IN£ compris dans ladite zone Zk>n, la distribution DDkn étant calculée le long d’une direction normale à la direction d’analyse Dk de ladite zone Zkn; pour chaque direction d’analyse Dk, le calcul d’un critère de variabilité CVk n de distributions {pDk n-v,DDk n_p+1,... ,DDk n_1,DDkn\ calculées pour des zones Zkin 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 n est optimum.
5. 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 {ZktTl} pour chaque image INn du flux d’images numériques {INn} , ledit calcul comprenant : le calcul d’un ensemble de rectangles 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 {Zk>n} comme étant égales respectivement à l’intersection des rectangles {RkiTl} avec l’image INn.
6. 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.
7. 7. Procédé selon la revendication 6, caractérisé en ce le petit côté de chaque rectangle Rkn de l’ensemble de rectangles {/?fe n} est au moins égal au diamètre du halo du disque solaire.
8. 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.
9. 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 (f?i,0£,Azr) du dispositif d’imagerie, définie par la position (θ[, Θ£) du capteur numérique par rapport à un plan horizontal et par la position d’un axe du capteur .par rapport à un azimut (Azr), et en ce que le procédé comprend une étape de détermination de la position prédéterminée θ^,Αζ7") du dispositif d’imagerie, comportant : l’acquisition d’un flux d’images numériques de calibration horodatées {/1V^} comprenant chacune le disque solaire ; pour chaque image 1N^ dudit flux {IN£} , le calcul d’une position P£ du centre du disque solaire dans l’image ΙΝξ en analysant l’image ΙΝξ ; le calcul de la position du dispositif d’imagerie {ΘΙ,Θ^,Αζ^ comme solution d’un problème d’optimisation selon la relation :

   où Ω est un ensemble de recherche, et || || est une fonction de coût quantifiant un écart entre deux vecteurs de M3.
10. 10. Procédé selon l’une quelconque des revendications 4 à 9, caractérisé en ce que le critère de de variabilité CVkn pour la direction Dk est calculé en fonction des distributions {DDkn_p,DDkn-p+1,....DD^^.DDkn] associées aux p> 2 dernières images numériques acquises.
11. 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 :

    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 nlème image et la (n-p),ème image dans le flux d’images numériques. {INn} , avec p strictement inférieur à n, £ référençant la klème direction Dk, et ||l/|| est une norme du vecteur V de l’espace W, /étant la dimension des distributions DDkn.
12. 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 DDkn, ledit masque est appliqué aux images du flux d’images numériques {INn} et/ou aux images segmentées {IN£} .
13. 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 fonction de ladite vitesse moyenne Vn et de ladite direction Dk n.
14. 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.
15. 15. Système d’estimation de la vitesse moyenne Vn de déplacement d’objets en mouvement selon une direction Dkoptin dans un flux d’images numériques {lNn} 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 scène; une unité de traitement d’information configurée pour : o l’estimation d’une direction Dkoptin en fonction du flux d’images numériques {INn} ; o pour chaque image INn du flux d’images numériques {INn} : le calcul d’une image segmentée IN% identifiant les objets en mouvement dans l’image INn; le calcul d’une distribution DVk 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 Dk n déterminée ; o le calcul de l’inclinaison d’au moins une structure affine comprise dans une matrice DVn = [DVkopt>n_qDVkopt>n_q+1... DV^^DV^^] issue de la concaténation des distributions DVkopt n des images entre la (n-q)ieme image du flux jusqu’à la nieme du flux d’images numériques {INn} , où q est un entier supérieur ou égal à 1 ; o et le calcul de la vitesse moyenne Vn en fonction de l’inclinaison calculée.
16. 16. Système selon la revendication 15 configuré pour mettre en œuvre un procédé selon l’une qpelconque des revendications 2 à 14.