FR 2998363 A1 20140523 FR 1260963 A 20121119 Procédé de détermination, dans un référentiel fixe 3D, de la localisation d'un engin en mouvement, dispositif et programme d'ordinateur associés La présente invention concerne le domaine de la détermination, dans un référentiel fixe 3D, de la localisation d'un engin mobile doté d'un capteur d'images associé à un premier référentiel local 3D, correspondant à la prise d'une image, par ledit capteur, de l'environnement de l'engin mobile comprenant au moins deux points caractéristiques dudit environnement. Certains de ces engins volants, comme les drones par exemple, sont utilisés pour accomplir des tâches d'exploration, reconnaissance, surveillance etc. Pour mener à bien l'ensemble de ces tâches, un tel engin volant doit être capable de connaître sa localisation, notamment pour réaliser des vols stationnaires. Dans de nombreux cas, il n'est pas possible pour l'engin volant de recourir à une localisation de type GPS (« Global Positioning system »). Il est donc nécessaire de le doter d'un outil lui permettant de réaliser un vol stationnaire sans GPS, en utilisant uniquement des mesures issues de capteurs inertiels et de capteur(s) d'images. Toutefois, toutes les approches proposées sont à base de filtres, ce qui signifie que des données initiales de localisation doivent être fournies. Or une erreur sur ces données d'initialisation peut corrompre de façon très importante l'ensemble des localisations estimées. En outre, ces approches nécessitent un grand nombre de points caractéristiques. Un point caractéristique est un point de l'environnement physique de l'engin volant, qui sert de repère dans les images, et qui est fixe (ou alors dont le mouvement est connu). A partir d'une image prise par le capteur d'images et comprenant un point caractéristique, la direction de ce point caractéristique est déterminée dans un repère tridimensionnel local à l'engin volant. L'utilisation d'un grand nombre de points caractéristiques génère une charge importante de calculs. L'une de ces approches connues utilise l'algorithme PTAM (en anglais « Parallel Tracking And Mapping ») et un filtre de Kalman étendu (cf. par exemple « Vision based navigation for micro helicopters », de Weiss, Stephan M., http://e- collection.library.ethz.ch/view/eth :5889). La présente invention vise à proposer une solution qui contribue à résoudre ces inconvénients de l'art antérieur. En effet, de façon non attendue, la Demanderesse a découvert qu'une fois que les valeurs de tangage et de roulis correspondant à une image ont été déterminées dans un référentiel absolu (i.e. dont la position est non dépendante du mouvement de l'engin volant), il est possible de déterminer la position de l'engin volant et la valeur de lacet correspondant à cette image, dans le référentiel absolu, et ce sans prendre en compte des valeurs d'initialisation de position ou de lacet. A cet effet, suivant un premier aspect, l'invention propose un procédé du type précité caractérisé en ce que les angles de roulis et de tangage de l'engin mobile dans le référentiel fixe correspondant à ladite prise d'image ayant été déterminés dans une étape préalable, les étapes suivantes sont mises en oeuvre, relativement à ladite image prise, pour déterminer l'angle de lacet et/ou les coordonnées de l'engin mobile dans le référentiel fixe : a/ on définit un deuxième référentiel local à l'engin mobile comprenant un axe présentant un écart d'angle connu avec l'axe de lacet du référentiel fixe, ledit deuxième référentiel local résultant de la succession de deux rotations, une desdites rotations autour d'un premier axe du premier référentiel local étant fonction de l'angle de roulis déterminé et l'autre desdites rotations étant fonction de l'angle de tangage déterminé ; b/ les coordonnées de chacun des deux points caractéristiques dans ledit deuxième référentiel local étant égales, à un facteur commun inconnu près, au résultat de la multiplication d'une matrice définissant ladite succession de rotations, par une matrice de paramètres intrinsèques du capteur d'images et par les coordonnées d'un pixel correspondant audit point caractéristique sur l'image, on détermine les valeurs prises par les quotients entre lesdites coordonnées de manière à s'affranchir du facteur commun inconnu ; c/ on exprime ensuite lesdites coordonnées des deux points caractéristiques dans ledit deuxième référentiel local en fonction de la distance entre lesdits points caractéristiques et desdites valeurs de quotients déterminées. Dans des modes de réalisation, le procédé suivant l'invention comporte en outre une ou plusieurs des caractéristiques suivantes : - on calcule les coordonnées, dans le deuxième référentiel local, du vecteur joignant les deux points caractéristiques, et selon lequel l'angle de lacet de l'engin mobile dans le référentiel fixe correspondant à ladite image prise est déterminé en fonction de la fonction arctan appliquée à une variable fonction desdites coordonnées, dans le repère local, du vecteur joignant les deux points caractéristiques ; - on calcule les coordonnées de l'engin mobile dans le référentiel fixe correspondant à ladite image prise, en fonction de coordonnées d'au moins un desdits points caractéristiques et dudit angle de lacet déterminé dans le référentiel fixe ; - étant donné trois points caractéristiques dans le plan perpendiculaire à l'axe de lacet du référentiel fixe, les étapes suivantes sont mises en oeuvre relativement à ladite image prise : - pour chacune des trois paires de points caractéristiques, on met en oeuvre les étapes b/ et cl et on détermine ainsi trois valeurs d'angles de lacet en fonction des valeurs de roulis et tangage déterminées dans l'étape préalable ; - à partir des trois valeurs d'angles de lacet déterminées, on met en oeuvre une méthode de minimisation d'une fonction cost() - (YW13 17W12 71)2 + (YW23 I7W12 72)2 dans laquelle les variables YW13, YW12 et YW23 sont des fonctions non-linéaires de = LRoulis,Tangage]T où 71 et 72 représentent deux angles internes du triangle de sommets desdits points caractéristiques ; - en fonction des valeurs de roulis et tangage issues de cette méthode de minimisation, les étapes b/ et cl sont mises en oeuvre relativement à deux des points caractéristiques et on détermine ainsi un angle de lacet dans le référentiel fixe. - le procédé comprend une étape de détermination des angles yi et 72 , selon laquelle on calcule : a c 2 b2 2 a c b 2 2 2 = arccos( ) et 72 - arccos( 2ab 2ac où a = V(/-12 )2 + (y2 - Y1)2 b = -à1,13 111)2 + (y3 - Y1)2 C = V(113 -112 )2 + (y3 - y2 )2 5 xi Yi les valeurs du/ =- et y =71, où z11 Yi sont les coordonnées du point Z- caractéristique P, dans le deuxième référentiel local, = 1 à 3, étant déterminées à partir xi = Rot.A-1. ui vi 1 des égalités suivantes : 2. , où u,, v, sont les coordonnées d'un pixel Yi zi_ de la représentation du point caractéristique P, dans l'image, A est une matrice prédéterminée de paramètres intrinsèques au capteur d'images, 2. un facteur de valeur inconnue ; et Rot est la matrice représentant la succession des deux rotations. - ladite étape de détermination des angles yi et 72 est mise en oeuvre dans une phase d'observation d'une méthode de filtre de Kalman ; - le référentiel fixe est défini tel que le premier point caractéristique en est l'origine, que l'axe de lacet est de direction contraire au vecteur gravité et que le deuxième point caractéristique est de coordonnée égale à zéro sur l'axe de tangage du référentiel fixe ; - un axe du deuxième référentiel local est choisi parallèle à l'axe de lacet du référentiel fixe ; - lorsque les deux points caractéristiques ont une même valeur de coordonnée sur l'axe de lacet du référentiel fixe, la valeur de coordonnée z de la localisation de l'engin mobile sur l'axe de lacet du référentiel fixe est : z=+ D avec A.,u = pi -#2;Ay=y1 -y2 et D la distance entre les deux V42 ±Av2 points caractéristiques, X1 V1 Y1 X2 Xi et X2 sont les coordonnées ill = - I V 1 = 7 /12 = - et4. 311 Z1 Y2 Z2 V2 = y2ou Z1 Z1 Z2 Z2 des deux points caractéristiques dans le deuxième référentiel local ; - lorsque les deux points caractéristiques n'ont pas une même valeur de coordonnée sur l'axe de lacet du référentiel fixe, la valeur de coordonnée z de la localisation de l'engin mobile sur l'axe de lacet du référentiel fixe est une solution de : (42 ± Av2)( z )2 2(4/12 ± Avv2) 42 z ± (A22 avec D D D) 1= 0 D A.,u = pi -,u2;4y = y1 -y2 et D la distance entre les deux points caractéristiques, X1 , ill X2 Xi et X2 sont les coordonnées /11 = I V 1 = 7 /12 = - et4. 311 Z1 Y2 Z2 V2 = ou Z1 Z1 Z2 Z2 des deux points caractéristiques dans le deuxième référentiel local. - le procédé comprend une étape de détermination de l'angle de lacet YW de l'engin mobile dans le référentiel fixe à l'aide de la formule suivante : YW = -0 = - 2 .arctan , \hz(/-12 - /11 ))2 + (- z(v2 - Vi) Z(/2 - y1 ))2 + (- ZU-12 - /11 b - le procédé comprend au moins une des étapes i/ et ii/ ci-dessous, avec -0 = - 2 .arctan - z(v2 - v1) .\/(4/-12 - /-11 »2 ± (Z(/2 - 1/1 ))2 ± (- ZU-12 - /-11 b i/ une étape de détermination de la coordonnée x de la localisation de l'engin mobile sur l'axe de roulis du deuxième référentiel fixe (Rn) en calculant x = cos e. z.,ui + sin e. z.v1 ; ii/ une étape de détermination de la coordonnée y de l'engin mobile sur l'axe de tangage (Y) du référentiel fixe en calculant y = -sin e. z.,ui + cos e. z.v1. Xi Y1 2 et _ Y2 - /-11 = 71, vi = , #2 _ x -72 et y2 - - Z2 sont calculées sur la base des égalités z1 suivantes : .1. Xi = Rot.A-1. u1 vi , où u1, vl sont les coordonnées d'un pixel de la 311 z1 1 représentation d'un des deux points caractéristiques dans l'image ; .1. x2 = Rot.A-1. u2 , où u2, v2 sont les coordonnées d'un pixel de la Y2 Z2 V2 1 représentation de l'autre des deux points caractéristiques dans l'image, A est une matrice prédéterminée de paramètres intrinsèques au capteur d'images 2. un facteur de valeur inconnue ; et Rot est la matrice représentant la succession des deux rotations ; Xi y1Xn Y2 - /11 = 71, vi =71, 112 = et V2 = sont calculées sur la base des égalités z2 Z2 suivantes : .1. Xi = Rot.A-1. u1 vi , où u1, vl sont les coordonnées d'un pixel de la 311 z1 1 représentation d'un des deux points caractéristiques dans l'image ; .1. x2 = Rot.A-1 u2 V2 , où u2, v2 sont les coordonnées d'un pixel de la Y2 Z2 1 représentation de l'autre des deux points caractéristiques dans l'image, A est une matrice prédéterminée de paramètres intrinsèques au capteur d'images (7), .1. un facteur de valeur inconnue ; et Rot est la matrice représentant la succession des deux rotations. Suivant un deuxième aspect, la présente invention propose un programme d'ordinateur pour localiser, dans un référentiel fixe 3D, un engin mobile doté d'un capteur d'images associé à un premier référentiel local 3D correspondant à la prise d'une image, par ledit capteur, de l'environnement de l'engin mobile comprenant au moins deux points caractéristiques dudit environnement, ledit programme comportant des instructions pour mettre en oeuvre les étapes d'un procédé suivant le premier aspect de l'invention lors d'une exécution du programme par des moyens de traitement, les angles de roulis et de tangage de l'engin mobile dans le référentiel fixe correspondant à ladite prise d'image ayant été déterminés dans une étape préalable. Suivant un troisième aspect, la présente invention propose un dispositif de détermination, dans un référentiel fixe 3D, de la localisation d'un engin mobile doté d'un capteur d'images associé à un premier référentiel local 3D correspondant à la prise d'une image, par ledit capteur, de l'environnement de l'engin mobile comprenant au moins deux points caractéristiques dudit environnement, ledit dispositif de détermination étant adapté pour déterminer préalablement les angles de roulis et de tangage de l'engin mobile dans le référentiel fixe correspondant à ladite prise d'image ; ledit dispositif de détermination étant caractérisé en que qu'il est adapté pour déterminer en outre l'angle de lacet et/ou les coordonnées de l'engin mobile dans le référentiel fixe en définissant un deuxième référentiel local à l'engin mobile comprenant un axe présentant un écart d'angle connu avec l'axe de lacet du référentiel fixe, ledit deuxième référentiel local résultant de la succession de deux rotations, une desdites rotations autour d'un premier axe du premier référentiel local étant fonction de l'angle de roulis déterminé et l'autre desdites rotations étant fonction de l'angle de tangage déterminé ; les coordonnées de chacun des deux points caractéristiques dans ledit deuxième référentiel local étant égales, à un facteur commun inconnu près, au résultat de la multiplication d'une matrice définissant ladite succession de rotations, par une matrice de paramètres intrinsèques du capteur d'images et par les coordonnées d'un pixel correspondant audit point caractéristique sur l'image, ledit dispositif de détermination étant adapté pour déterminer les valeurs prises par les quotients entre lesdites coordonnées de manière à s'affranchir du facteur commun inconnu ; ledit dispositif de détermination étant adapté pour exprimer ensuite lesdites coordonnées des deux points caractéristiques dans ledit deuxième référentiel local en fonction de la distance entre lesdits points caractéristiques et desdites valeurs de quotients déterminées. Ces caractéristiques et avantages de l'invention apparaîtront à la lecture de la description qui va suivre, donnée uniquement à titre d'exemple, et faite en référence aux dessins annexés, sur lesquels : - la figure 1 représente un engin volant ; la figure 2 représente un dispositif de traitement dans un mode de réalisation de l'invention ; la figure 3 est un organigramme de différentes étapes mises en oeuvre dans un mode de réalisation de l'invention ; la figure 4 représente des points dans des repères tridimensionnels utilisés dans un mode de réalisation de l'invention ; la figure 5 est une vue de deux points caractéristiques dans un plan (X, Vn) dans un mode de réalisation de l'invention ; la figure 6 est un organigramme de différentes étapes mises en oeuvre dans un mode de réalisation de l'invention ; la figure 7 représente des points dans des repères tridimensionnels utilisés dans un mode de réalisation de l'invention. En référence à la figure 1, on caractérise généralement la localisation d'un engin volant 10 d'une part par la détermination de la position, dans l'espace, de son centre de gravité O, et d'autre part par son attitude, c'est-à-dire par ses rotations, mesurées à l'aide d'angles de rotation, par rapport à trois axes définissant un repère tridimensionnel passant par un point d'origine du repère, par exemple le centre de gravité O : l'angle de roulis (en anglais « roll ») par rapport à un axe de roulis AR, l'angle de tangage (en anglais « pitch ») par rapport à un axe de tangage AR et l'angle de lacet (en anglais « yaw ») par rapport à un axe de lacet Ay. Dans le cas représenté en figure 1, l'axe de roulis AR considéré est l'axe longitudinal de l'engin volant 10, l'axe de tangage AR est l'axe transversal de l'engin volant et l'axe de lacet Ay est l'axe perpendiculaire aux axes AR et A. Dans le mode de réalisation considéré de l'invention, un engin volant 10, par exemple un drone, comporte un dispositif de traitement 1, en référence à la figure 2, adapté pour effectuer des mesures de données inertielles et des captures d'images et pour traiter ces données. Ce dispositif de traitement 1 comporte un microprocesseur 2, une mémoire de stockage de données 3, un ensemble 4 de capteurs inertiels et un capteur d'images 7. Dans le mode de réalisation considéré, l'ensemble 4 de capteurs inertiels comporte un module 5 d'accéléromètre avec 3 accéléromètres orthogonaux entre eux, fournissant toutes les périodes de longueur TA, les mesures d'accélération linéaire du drone 10 relativement à trois axes orthogonaux (X1, Y1, Z1) passant par un point OL du drone et définissant ainsi, en référence à la figure 4, un repère R1= (OL, X1, YI, Z1), local au drone, et nommé ci-après repère local. Dans le mode de réalisation considéré, l'ensemble 4 de capteurs inertiels comporte en outre un module 6 de gyroscope comportant trois gyroscopes orthogonaux entre eux, qui fournissent, toutes les périodes de longueur TB, des mesures de rotation du drone 10 autour des axes (X1, YI, Z1). Le capteur d'images 7, par exemple une caméra monoculaire, est fixé sur le drone de manière à capturer une image de l'environnement du drone 10 toutes les périodes de longueur T. Les images capturées de l'environnement comportent notamment un nombre N de points caractéristiques de l'environnement. Dans le cas présent, l'axe Z1 est l'axe optique du capteur d'image 7 et le point OL est le centre optique du capteur d'images 10 7. De manière connue (cf. par exemple le lien Internet http://opency.willowgarage.com/documentation/camera calibration and 3d reconstructio n.html), les coordonnés d'un point de l'espace peuvent être retrouvées, à un facteur d'échelle près, en fonction des coordonnées du pixel représentant ce point dans une image capturée par le capteur d'images 7, à l'aide d'une matrice A prédéterminée de paramètres intrinsèques à ce capteur d'images, et indépendants de l'image. Ainsi pour un point M ayant pour coordonnées (xi, yl, z1) dans le repère local RI et représenté par un pixel de coordonnées (u,v) sur une image prise par le capteur d'images 7, on a: u y = .1.A. xi , où 2 est un facteur de proportionnalité positif, indépendant de 1 y' zi l'image . On obtient donc les coordonnées d'un vecteur xp Y p , proportionnel à xi , donnant ZP_ Y1 zi la direction avec du point M depuis u l'origine X =2 Xi = A-1. V OL du repère RI : P X .11/ z' YP 1 (1) yp P X I et - YI Z = - (2) P z1 Z1 zp zP Un point caractéristique est par exemple un coin d'angle entre deux structures de l'environnement du drone 10, le centre d'une balle, le centre d'un spot laser etc. Dans un mode de réalisation, TA = TG = 0,01 seconde (s), et TI = 0,1 s.. Dans la mémoire 3 du dispositif de traitement 1, un programme P d'instructions logicielles est stocké, qui, lorsqu'il est exécuté par le microprocesseur 2, est adapté pour mettre en oeuvre les étapes indiquées ci-dessous en référence à la figure 3. Dans un mode de réalisation considéré ci-dessous, N est choisi au moins égal à 2. Deux points caractéristiques considérés sont nommés P1 et P2. Dans une première étape 100, en référence aux figures 3 et 4, et étant donné ces deux points caractéristiques Pl, P2, un repère de référence absolu RG est défini. Dans le mode de réalisation considéré, le point caractéristique P1 est l'origine de ce repère de référence RG, défini en outre par trois axes orthogonaux X, Y et Z passant par P1 et tels que l'axe Z est l'axe parallèle au vecteur de gravité g (tel que déterminé par le module 5 d'accéléromètre) et de direction opposée, et le point caractéristique P2 est dans le plan (X,Z), i.e. le point caractéristique P2 a une valeur de coordonnée nulle sur l'axe Y. La détermination d'au moins les 6 valeurs définissant, dans le repère de référence absolu RG, la position (x, y, z) du point OL du drone et les angles de roulis R (relativement à l'axe X), tangage P (relativement à l'axe Y), lacet YW (relativement à l'axe Z) du drone 10, est réalisée dans les étapes 101 et 102. Dans l'étape 101, ces valeurs de roulis R et de tangage P « absolues » sont déterminées, par rapport à l'axe de roulis pris égal à l'axe X et à l'axe de tangage pris égal à l'axe Y du repère de référence RG. Plusieurs solutions peuvent être mises en oeuvre pour cette étape 101. Deux d'entre elles seront explicitées plus loin. Deux points caractéristiques minimum sont nécessaires pour réaliser l'étape 102, i.e. pour déterminer l'angle de lacet et la position par rapport au repère absolu RG. Lorsque plus de deux points caractéristiques sont disponibles, le même nombre d'images est requis que dans le cas de deux points caractéristiques, mais la précision dans les valeurs d'angles estimées est supérieure. On considère à présent une image capturée par le capteur d'images 7 et on lui associe une valeur R de roulis et une valeur de tangage T du drone correspondant à l'instant de la prise de cette image (en prenant directement les valeurs déterminées à l'étape 101 pour cette image ou en moyennant plusieurs valeurs de roulis, respectivement tangage, déterminées pour plusieurs images prises dans une plage de temps donnée fonction de l'instant de la prise de cette image considérée etc.). Dans l'étape 102, la valeur de lacet YW et la (position (x,y,z) du drone 10 dans le repère absolu RG sont déterminées. Comme indiqué précédemment, l'image considérée comporte une représentation de chaque point caractéristique P1 et P2, ce qui permet de déterminer la position de ces points (moyennant un facteur d'échelle) dans le repère local RI (cf. notamment les équations (1) et (2)). La connaissance des valeurs absolues de roulis et tangage, R et P, permet de définir un nouveau repère local R, et d'exprimer la position des points caractéristiques dans ce nouveau repère local Rn. Ainsi, en référence à la figure 6, dans une étape 102 1, on définit ce nouveau repère local Rn = (OL, Xn, Yn, Zn), qui a pour origine le point OL et comporte 3 axes orthogonaux Xn, Yn, Zn. Il est représenté en pointillés sur la figure 4. L'axe Z, est défini de sorte qu'il a la même orientation que l'axe Z du repère général RG. L'axe n représente la projection de l'axe Z, sur le plan (YI, Z1). L'angle de tangage P est l'angle entre l'axe Z, et l'axe n, et l'angle de roulis R est l'angle entre l'axe n et l'axe Z1. Le but est de déterminer les coordonnées Y z du point OL dans le repère absolu RG, ainsi que l'angle YW de l'axe Xn par rapport à l'axe X (YW étant la valeur de l'angle de lacet dans le repère absolu RG). X2 On appelle y1 les coordonnées du point caractéristique P1, et y2 celles du zi z2 point caractéristique P2, dans le nouveau repère local R. On nomme pi = et v Y2 5 1_ 5 -//2_ 2 _ - Z1 Z1 Z2 Z2 A l'aide de l'équation (1), à partir de l'image considérée, et en fonction en outre d'une matrice de rotation, nommée Rot, définissant l'opération de rotation permettant à partir du repère local RI d'obtenir le nouveau repère local Rn, on a : 2. y1 = Rot.A-1. u1 , où u1, vl sont les coordonnées du pixel sur lequel la zi 1 représentation du point P1 dans l'image est centrée; et = Rot.A-1. , où u2, v2 sont les coordonnées du pixel sur lequel la X2 2. Y2 z2 U2 V2 1 représentation du point P2 dans l'image est centrée; et Rot = cos acos /3 cos a sin /3 sin y- sin acos y cos a sin /kas y+ sin a sin y sinacosfl sinasinfisin y+cosacos y sinasinficosy-cosasin y -sinfi cos fisin y cos /kas y Où: a est l'angle de rotation effectué autour de l'axe Z, correspondant au lacet et fixé à 0, fi est l'angle de rotation effectué autour de l'axe YI et correspondant au tangage , et y est l'angle de rotation effectué autour de l'axe X, correspondant au roulis, pour transformer le repère RI en le repère Rn. La matrice Rot est obtenue à l'aide d'une première rotation d'une valeur d'angle y autour de l'axe XI, suivie d'une deuxième rotation, de la valeur d'anglet autour de l'axe YI , et enfin d'une troisième rotation, de la valeur d'angle a autour de l'axe Z1 (donc (3 = Pet y = R calculés à l'étape 101). Pour obtenir le repère Rn, seules les deux premières rotations doivent être réalisées. En d'autres termes, les deux repères Rn et RI ont la même valeur de lacet (inconnue). Ceci revient à fixer la valeur de a à 0 dans la matrice Rot, (pour plus d'information sur cette matrice de rotation, cf. par exemple http://planning.cs.uiuc.edu/node102.html). En calculant Rot.A-1. u1 vi et Rot.A-1 u2 V2 , on en déduit les valeurs de /il = ), 1 1 z1 Y1, x2 ,f , Y2 V1 = 5 i-42 _ - - 2 _ - _ - Z2 '' Z1 2 Z2 On cherche à déterminer les coordonnées X du point OL dans le repère absolu Y z RG, ainsi que l'angle de lacet YW par rapport à l'axe Z de ce repère RG. Dans une étape 102 2, on détermine la valeur de z. Cas 1 : Pl P2 sont tous deux dans le plan (X, Y) Dans une étape 102 21, qui est mise en oeuvre lorsque les points caractéristiques appartiennent tous deux au plan (X,Y), i.e. ont la même coordonnée sur l'axe Z du repère absolu RG, on détermine la valeur de la coordonnée z du point 01_ sur l'axe X, avec z1 = z2= -z : Avec P1 = -z. Ili et P2 = -z 112 V2 et en appelant D la distance séparant P, de P2, VI 1 1 il s'ensuit que z=± D; Ail = iii -/-12;Av =v1- v2 . VA,I2 ± Av2 En d'autres termes, z est exprimée en fonction de D, au signe près. Le signe de z est ensuite déterminé. En effet, le signe de la coordonnée sur l'axe Z, des points caractéristiques Pl, P2 est positif (sinon ils ne seraient pas sur l'image). On en déduit le signe de la coordonnée (-z) des points caractéristiques Pl, P2 sur l'axe Z, ou Z (la matrice Rot dépendant des valeurs, connues, de roulis et tangage). On en déduit x1 = -z.,ui yi = -z.vi x2 = -z.,u2 et y2 = -z.v2 . En référence à la figure 5, dans une étape 102 3, on en déduit la valeur du lacet YW dans le repère absolu RG à l'aide de la fonction arctan2 : y' - YW = -0 = - arctan2(y2-yi,x2-x1) = - 2 .arctan I Y2(4) .N/ (x2 -x1)2 + (y2 -Y1)2 + ( ,x2 -.X1) De la figure 5, on déduit en outre les coordonnées x et y de l'origine OL des repères locaux RI et Rn dans le repère global RG : X = -COS 0. Xi - sin a y, (5) y = sin 0. xi - cos a y, (6) X Les coordonnées Y z de 01_ représentant la position du drone dans le repère absolu RG sont ainsi déterminées, ainsi que l'angle de lacet YW recherché par rapport à ce repère absolu. Elles sont exprimées en fonction de la distance D. Les coordonnées de cette position dépendent proportionnellement de D. Cas 2 : P2 n'est pas dans le plan (X, Y) Alors z2 # z1 et z1 = -z et z2 = -z + Az, avec Az #0. Dans une étape 102 22, qui est mise en oeuvre lorsque les points caractéristiques n'appartiennent pas tous deux au même plan (X,Y), i.e. ont la même coordonnée sur l'axe Z du repère général RG, on exprime les coordonnées de P1 et P2 de la façon suivante cette fois en vue de déterminer la valeur de la coordonnée z du point 01_ sur l'axe X : Ili 112 P1 = -z. VI et P2 = (-z + Az ). V2 et D étant la distance séparant P1 de P2 et 1 1 Adu =,u1-,u2;Ay = v1 - v2 , on obtient l'équation polynomiale du second ordre en z : (A,u2 ± A v 2)(-Z)2 - 2(41/2 ± A vv2)Az-z +(-Az)2 -1 = 0 D D D D Si le rapport D Aest connu, l'équation précédente peut être utilisée pour déterminer -z . Elle fournit deux solutions. De la même façon que dans le cas 1, une de ces solutions D est éliminée du fait que le signe de zl, z2 est positif. En d'autres termes, z est exprimée en fonction de D, au signe près. Az On en déduit x1 = -z.iii , y1 = -z.vi , x2 = z(-1+ -).1u2 et y2 = -z(-1+ Az).v2. D D Puis, l'étape 102_3 est mise en oeuvre comme dans le cas 1 et on déduit la valeur du lacet YW dans le repère absolu RG à l'aide de la fonction arctan2 : YW = -0 = - arctan2(y2-Y,,x2-x1) = - 2 .arctan V(x2 -x1)2 + (y2 -y1)2 + (x2 -x1) De la figure 5, on déduit en outre les coordonnées x et y de OL dans le repère absolu RG : X = -cos 0. x1 - sin 0. yl (5) y= sin 0. x,- cos a yl (6) Y2 Yi (4) Les coordonnées x de 01_ représentant la position du drone dans le repère Y z absolu RG sont ainsi déterminées, ainsi que l'angle de lacet par rapport à ce repère absolu. Elles sont exprimées en fonction de la distance D. Les coordonnées de cette position dépendent proportionnellement de D. Ainsi, selon l'invention, les valeurs de position (x, y, z) et de lacet du drone sont déterminables selon l'invention dans un repère absolu à l'aide de deux points caractéristiques seulement et d'une image pour laquelle les valeurs de roulis et de tangage ont préalablement été déterminées, sans besoin de valeurs de position d'initialisation ou de valeur de lacet d'initialisation, et avec un nombre réduit de calculs, permettant des missions autonomes du drone. Dans une étape 103, facultative, une fois la position du drone et les angles de tangage, roulis et lacet déterminés dans le repère absolu RG, des traitements et/ou actions sont réalisées en fonction de ces éléments. Dans un mode de réalisation, la valeur de D peut être estimée, par exemple en intégrant des données d'images fournies par le capteur d'images 7 et en fonction des données inertielles fournies par le module d'accéléromètre 5 et le module de gyroscope 6, obtenues lorsque le drone se déplace avec une accélération non nulle, avant de passer en vol stationnaire (cf. notamment le document mentionné plus haut A. Martinelli, Vision and IMU Data Fusion: Closed-Form Solutions for Attitude, Speed, Absolute Scale and Bias Determination, Transaction on Robotics, Volume 28 (2012), Issue 1 (February), pp 44-60). On notera que même si la valeur de D est obtenue avec une erreur non négligeable (par exemple de 10%), cette erreur n'affectera pas la performance de vol stationnaire. En effet, la solution décrite fournit la position comme fonction linéaire de D. Donc l'erreur sur D affectera de façon systématique la position. Ce serait différent si la solution déterminait la position avec une erreur statistique d'amplitude 10%. Il serait alors impossible de stabiliser le drone puisqu'à chaque pas de calcul, la position courante estimée différerait de façon significative de la précédente. Selon la solution décrite, la différence entre des estimations de position successives est indépendante de l'erreur sur l'estimation de D. Une erreur dans les estimations des valeurs de roulis et de tangage effectuées à l'étape 101 affecte de façon linéaire les estimations de la position et du lacet telles que décrites ci-dessus en référence à l'étape 102. Dans un mode de réalisation, pour réduire l'impact d'une erreur d'estimation des valeurs de roulis et de tangage sur les estimations de la position et du lacet, l'étape 102 est remplacée par une étape 102bis, dans laquelle trois points caractéristiques P1, P2 et P3 au moins apparaissant sur l'image considérée et faisant partie du plan (X,Y), sont pris en compte, au lieu des deux points caractéristiques P1, P2 considérés dans l'étape 102. Sur la figure 7, les trois points caractéristiques Pl, P2 et P3 sont représentés. Ces points Pl, P2 et P3 sont coplanaires, dans le plan (X,Y) du repère RG. Ils forment un triangle. On considère que les angles internes de ce triangle sont connus, ce qui correspond à la connaissance des angles yi et 72 représentés en figure 7. Si ce n'est pas le cas, une étape d'estimation préalable de ces angles devra être réalisée, comme décrit plus loin. Dans l'étape 102bis, ces opérations décrites ci-dessus dans l'étape 102 sont mises en oeuvre relativement aux référentiels RG, Rn, RI, en considérant dans un premier temps les deux points Pl, P2, tous deux dans le plan (X,Y), dans un deuxième temps en considérant les deux points Pl, P3, tous deux dans le plan (X,Y), et dans un troisième temps en considérant les deux points P2, P3, tous deux dans le plan (X,Y), comme détaillé ci-dessous xi Yi les coordonnées du point caractéristique Pi, i = 1 à 3, dans le nouveau zi_ repère local R, et du/ = vi = zi zi En mettant en oeuvre l'étape 102 dans un premier temps en considérant les points caractéristiques Pl, P2, on obtient la valeur de la coordonnée z du point OL et une première estimation de lacet YW12 (correspondant à la valeur de la rotation du drone autour de l'axe Z par rapport à l'axe X = (P1P2)) `IM 12 = - arctan2(y2-y, ,x2-x1) D et z=± ; en appelant D la distance séparant P, de P2. 11(/-11 - /12 )2 + (Y1 - y2 )2 Comme vu dans la description de l'étape 102, le signe de z et les valeurs de du, et v, , i= 1 à 3 sont fournis par les données du capteur d'images 7. En mettant en oeuvre l'étape 102 dans un deuxième temps en considérant les points caractéristiques Pl, P3, on obtient une deuxième estimation de lacet YW13 pour l'image considérée (correspondant à la valeur de la rotation du drone autour de l'axe Z par rapport à l'axe (P, P3)) YW13 = arCtan2(Y3-Y1,X3-)(1). En mettant en oeuvre l'étape 102 dans un troisième temps en considérant les points caractéristiques P2, P3, on obtient une troisième estimation de lacet YW23 pour Soit l'image considérée (correspondant à la valeur de la rotation du drone autour de l'axe Z par rapport à l'axe (P2P3)) YW23= - arctan2(y3-y2,x3-x2)- Comme la figure 7 le montre, les trois angles de lacet calculés YW12, YW13, YW23 sont en outre liés par les équations suivantes : YVV13= YW12+ YW23= YW12+ 72 - On en déduit que : = YW13 YW12 72 = YW23 -YW12 - Comme les angles yi et 72 sont connus, par exemple yi = 2;1- 3 et Ï2 = 3 , l'estimation de l'angle de lacet recherchée, ie YW=YW12, peut être corrigée, en prenant en compte cette contrainte. Cette prise en compte peut être effectuée en définissant une fonction de coût, nommée cost, et en cherchant alors à minimiser la fonction de coût, par exemple à l'aide de la méthode de Levenberg-Marquardt. Par exemple, la fonction de coût est définie à l'aide d'une méthode des moindres carrés non linéaire. Dans un mode de réalisation, la fonction cost, définie comme la différence quadratique entre les valeurs attendues (dépendant des paramètres modélisés) et les valeurs connues est définie ainsi : cos t() = (Y1/4/13 - Y1/4/12 - Y1 )2 + ( YW23 - YW12 - 72)2 , dans laquelle les variables YW,, sont des fonctions non-linéaires de [Roulis,Tangage]T ]. Une fois que les angles de Roulis et Tangage qui minimisent la fonction de coût cost sont déterminés, qui sont donc les angles de roulis et tangage corrigés, l'angle de lacet corrigé YW, puis les coordonnées x y, et z de OL sont calculés en fonction de cet angle de lacet corrigé YW, à partir des équations (5) et (6) détaillées ci-dessus en référence à l'étape 102. Si les angles internes du triangle (P1, P2, P3) ne sont pas connus, une étape d'estimation préalable des angles yi et 72 doit être mise en oeuvre, par exemple à l'aide d'un filtre de Kalman comme décrit ci-dessous. L'état à estimer est F = (71, 72 )7- - Des valeurs d'initialisation sont choisies pour les angles yi et 72 (par exemple celles indiquées plus haut yi =3 3 et 72 = ) , et pour la matrice de covariance initiale -2 ir Pt,,, (par exemple, la valeur initiale pour Ptri est 6p2t 1 0 , avec la valeur (Ypt étant égale 0 1 à l'incertitude en radians sur la connaissance des valeurs d'angle yi et 72 . Pendant la phase de prédiction, le filtre de Kalman ne met pas à jour l'état F, ni la matrice de covariance, car les angles yi et 72 sont constants dans le temps. Pendant la phase d'observation, les valeurs de roulis et tangage estimées à l'étape 101 et les images, sont utilisées. Connaissant les valeurs de roulis et tangage estimées, une rotation du repère local RI en le nouveau repère local R, peut être réalisée, avec une erreur due à l'estimation de ces angles de roulis et tangage, à l'aide de la matrice Rot. x, Les coordonnées Y des points caractéristique Pi, i = 1 à 3 dans le nouveau Z - repère local R, sont telles que : x Y =z vi z 1 Les données d'image et les matrices Rot et A fournissent les valeurs de du, et v, , comme vu précédemment, i. = 1 à 3. On peut donc ensuite calculer les longueurs a, b, c des côtés du triangle (P1, P2, P3) où a = II PIPA, b = °PIPA et c = II P2P3 II: a = z'(u2 - )2 + (1/2 - )2 b = 4(113 - 111)2 ± (v3 - v 1)2 C = 4(1/3 -/2 )2 + (v3 v2 )2 Par application de la loi du cosinus : c2 a2 ^ b2 2 2abcos yi et b2 _ a2 -2 ^ C 2accos(ir - y2 ) . La valeur des angles yi et 72 peut ainsi être calculée à l'aide des formules 2 suivantes : = arccos(a +b c2 ) et 72 = arccos( a2 + 2 b2 2ab 2ac On notera que les valeurs de ces angles yi et 72 sont calculées de façon indépendante de z. On peut alors utiliser les valeurs calculées de yi et 72 comme observation directe de l'état F du filtre de Kalman pour corriger la valeur précédente de F. L'estimation des valeurs des angles converge ainsi vers les valeurs effectives yi et 72 - Dans les modes de réalisation décrits ci-dessus, le repère local considéré est celui du capteur d'images 7. D'autres repères locaux peuvent être sélectionnés, l'équation (1) donnant la correspondance entre les coordonnées d'un pixel sur une image et les coordonnées dans le repère local devant alors être adaptées en conséquence de façon connue (cf. http://opency.willowgarage.com/documentation/camera calibration and 3d reconstructio n.html) L'invention a été décrite en référence à un drone ci-dessus, mais elle peut bien sûr être mise en oeuvre pour tout type d'engin dont la localisation doit être déterminée, notamment des engins volants. Par ailleurs, comme indiqué plus haut, l'étape 101 peut être mise en oeuvre de diverses façons. Par exemple, une sous-étape d'estimation du vecteur de gravité dans le repère local RI est mise en oeuvre sin P Gx On nomme G le vecteur de gravité, G = g - sin Rcos P - cos Rcos P Gy G z dans le repère RI, où g - 9,8 ms-2 est l'amplitude de la gravité, P l'angle de tangage et R l'angle de roulis. Lorsque les composantes Gx, Gy et G, du vecteur G sont connues, on en déduit les valeurs P et R : P = arcsinGx et R = arctan 2(Gy,Gz ) si cos P 0. Ainsi, l'étape d'estimation des valeurs de roulis R et tangage P se ramène ici à une estimation de la gravité dans le repère local RI. On distingue deux cas : un cas où l'accélération inertielle est négligeable par rapport à l'accélération gravitationnelle, et un cas où l'accélération inertielle n'est pas négligeable par rapport à l'accélération gravitationnelle. Dans le cas où l'accélération inertielle est négligeable par rapport à l'accélération gravitationnelle, ce qui est fréquemment le cas, les composantes retenues pour le vecteur G sont celles fournies directement par le module 5 d'accéléromètre. Elles sont en outre éventuellement filtrées à l'aide d'un filtre de Kalman étendu (EKF). Dans ce filtrage, dans une phase de prédiction, le vecteur de gravité estimé subit une rotation en utilisant les données du module de gyroscope 6, puis une phase d'observation est mise en oeuvre en utilisant l'accélération fournie par le module d'accéléromètre 5, qui coïncide avec la gravité puisque l'accélération inertielle est supposée négligeable. Dans le cas où l'accélération inertielle n'est pas négligeable par rapport à l'accélération gravitationnelle, les composantes retenues pour le vecteur G sont calculées comme détaillé dans le document A. Martinelli, « Vision and IMU Data Fusion: Closed- Form Solutions for Attitude, Speed, Absolute Scale and Bias Determination », Transaction on Robotics, Volume 28 (2012), Issue 1 (February), pp 44-60 et rappelé ci-dessous. Cette méthode permet de déterminer plusieurs quantités physiques en collectant des données délivrées par l'ensemble de capteurs inertiels 4 et par le capteur d'images 7 pendant un court intervalle de temps Tful et permet d'exprimer la valeur de ces quantités physiques au temps seulement en fonction des données inertielles et des données d'images obtenues pendant l'intervalle [Tin, TA. Parmi ces quantités physiques figurent les trois composantes du vecteur gravité dans le repère local RI, i.e. les composantes Gx, Gy et G, de G. On considère le repère absolu RG défini ci-dessus à l'étape 100 (on notera toutefois que tout repère (x, y, z) dont l'axe z pointe verticalement vers le haut (i.e. est parallèle, et est opposé, au vecteur gravité) peut être considéré pour déterminer les valeurs de roulis et tangage dans le mode de mise en oeuvre ici décrit de l'étape 101. On utilisera ci-après des lettres minuscules en gras pour identifier les vecteurs dans ce repère RG (par exemple, le vecteur gravité sera repéré par g, avec g = où g est la norme du vecteur gravité). On suppose ici que le repère de l'ensemble de capteurs inertiels 4 et celui du capteur d'images 7 coïncident entre eux, et forme le repère dit local RI. On utilisera ci-après des lettres majuscules en gras pour identifier les vecteurs dans ce repère RI. Comme le repère local RI est dépendant du temps, on adopte la notation suivante : le vecteur Wt( ) dont les coordonnées varient en fonction du temps représenté par la variable z est, dans le repère local dans son état considéré à un instant t, le vecteur de coordonnées absolues w( ). En outre, on note C 12 la matrice qui caractérise la rotation du repère local RI survenue dans l'intervalle de temps [tt, t2]. On a alors : Wt, ( )= C W t 2( ) et ( Ct? )-1= C» Enfin Ct représente la matrice de rotation entre le repère absolu RG et le repère local RI au temps t, i.e. W(2 )=Ct W t(2 L'ensemble de capteurs inertiels 4 fournit des données d'accélération, i.e le vecteur d'accélération A, et des données de vitesse angulaire du drone 10. Le module d'accéléromètre 5 mesure non seulement l'accélération inertielle, i.e. le vecteur A'"rt'elle mais aussi l'accélération gravitationnelle, i.e. le vecteur G. On suppose que les images du capteur d'images 7 comprennent les représentations d'un point caractéristique, par exemple Pl, ou de plusieurs points caractéristiques pendant l'intervalle [T,,, Tful comme représenté en figure 4. La position r du drone à un instant t c [Tin, 'fin] satisfait l'équation suivante : t r(t)= r(Tin) + v(Tin)At + f f ag)gclz- (10) T ' T ' où v adr/dt, a a dv/dt et 31 - T. Le module d'accéléromètre 5 ne fournit pas le vecteur d'accélération a( mais il fournit l'accélération dans le repère local qui comporte aussi la composante d'accélération gravitationnelle. Plus particulièrement, le module d'accéléromètre 5 fournit Ad) (Ç) - On notera que la gravité est retranchée car, lorsque le drone 10 n'accélère pas (i.e. àetinertielleg) est nulle), le module d'accéléromètre 5 détecte une accélération qui est celle d'un objet accéléré vers le haut en l'absence de gravité. On notera aussi que le vecteur dépend du temps uniquement parce que le repère local peut tourner. L'équation (10) est réécrite à l'aide du vecteur Ae) fourni par le module d'accéléromètre 5 : r(t)= r(Tin) + v(Tin)At + g 22 + CTrn S, (t), t r où S (t)q. f CTn Ae)d (:1' " La matrice CT peut être obtenue en fonction de la vitesse angulaire du drone 10, pendant l'intervalle [Tin, telle que fournie par le module de gyroscope 6 (cf. J. A. Farrell, « Aided Navigation: GPS and High Rate Sensors », McGraw- Hill, 2008). Et le vecteur S (t) peut être obtenu par intégration des données fournies par le module de gyroscope 6 et le module d'accéléromètre 5 pendant l'intervalle [Tin, fl. Supposons que N points caractéristiques P,, i = 1 à N, soient observés simultanément. Selon la convention de notation retenue, la position du point caractéristique P, dans le repère absolu est représentée par le vecteur p', i = 1 à N et la position du point caractéristique P, au temps t, dans le repère local considéré au temps t, est représentée par le vecteur P't(t), i = 1 à N. Nous avons : = r(t) + CT Pt(t), et en utilisant l'équation (10), on obtient l'équation suivante, pour i = 1 à N : Pit(t) = CT" [P (Tin) )- V (Tin )4t- G 7- in S rn (t) ], Une image unique capturée au temps t fournit tous les vecteurs Pt(t) à un facteur d'échelle près. Comme les données fournies par le module de gyroscope 6 pendant l'intervalle [Tin, 'fin] permettent de construire la matrice CTn , obtenir les vecteurs Pt(t) à un facteur d'échelle près permet de déterminer aussi P (t) = C7-t Pt(t) à un facteur d'échelle près. De l'équation 11, on déduit, pour i = 1 à N : P = PIT', (Tin) - V Ln (Tin)A A2 t - G 7- n - 2t S (12) Supposons que la caméra fournisse ni images des mêmes N points caractéristiques à des instants consécutifs j = 1 à ni et t1 = Tin Pour simplifier la notation, à partir de ce point, la notation suivante est adoptée : P; = PT;n(1-1) ,i= 1 àN;j= 1 àni; P' (Tin) ,i=làN; V=VTn(Tin) G 'GT' ; Si = ST', (t ),j= 1 àni En outre, on nomme p ij le vecteur unité qui a la même direction que Pi. Et on introduit les inconnues À ij telles que P; = À . p Sans perte de généralité, on fixe Tin = 0, i.e.Ot = t. L'ensemble de capteurs inertiels 4 fournit p ij et Si pour i = 1 à N ; j = 1 à ni. L'équation (11) peut être écrite de la manière suivante : P' - Vt; - Gti2/2 - À tj .p j = (12) Le problème est la détermination des vecteurs P', i = 1 à N, V et G. On peut utiliser l'équation (12) pour déterminer ces vecteurs, ce qui nécessite de déterminer aussi les valeurs À'j . En considérant j= 1 dans l'équation (12), i.e. = t1 = = 0, on obtient facilement P' = .p . On peut alors écrire le système linéaire de l'équation (12) comme le système des deux équations (13) qui suit : _ GtJ2/2 - Vt; + À .p -À lj .p lj = À11 1 À1 " 1 À +À", .pj =03 1 " 1 j " I 1 où j = 2 à ni, i = 2 à N et 03 est le vecteur nul de dimensions 3 x 1. Ce système d'équations linéaires (13) consiste en 3(ni -1)N équations à N.ni + 6 inconnues. Définissons les vecteurs à deux colonnes X et S, et la matrice (1): X= [GT, VT, À , À iN Àln, ÀNn]1- 03r- 03 03 S =[S2T , 03,-. , 03, S3T , 03,-. , 03,-. , S ni T, 03,..., 03 03 03 T2 S2 /11 03 03 "1 r,3 03 03 03 033 033 1/11 112 03 03 Uni i 03 03 - - - 1 2 "1 "2 Pm; "3 " IV /12 1-12 03 I" ni 0 3 N 033 033 1/11 03 .1-41 "1 r, /-4"N 1 " 03 03 /-42 03 2 U ni kti 3 03 03 03 12 03 03 03 033 033 /11 /11 03 03 03 03 "1 033 033 1"1 r, 2 3 Où T , SJ 4,13 et 13 est la matrice identité de dimensions 3 x 3 ; 033 est la 2 matrice nulle de dimensions 3 x 3 (on note que le troisième ensemble de colonnes disparaît en cas d'absence d'erreur). Le système d'équations linéaires (13) s'écrit ainsi de façon compacte : cP X = S (14). Les données de l'ensemble de capteurs inertiels 4 sont contenues dans ce système linéaire. En outre, on suppose que l'amplitude de l'accélération gravitationnelle est a priori connue. Cette information supplémentaire est ajoutée dans le système linéaire via l'équation quadratique suivante : I G I = g. En introduisant la matrice FI de dimensions 3 x (N.n,+6), avec n.[I3, 03, ..., 03], cette contrainte quadratique peut être écrite en fonction de X de la manière qui suit : I nx 12 = g2. (15). Le problème est donc résolu en déterminant le vecteur X qui satisfait les équations (14) et (15). Dans le document mentionné ci-dessus, A. Martinelli, « Vision and IMU Data Fusion : Closed-Form Solutions for Attitude, Speed, Absolute Scale and Bias Determination », Transaction on Robotics, Volume 28 (2012), Issue 1 (February), pp 4460, il est montré qu'il existe une unique solution, deux solutions distinctes ou une infinité de solutions. Quand la matrice (1) est de plein rang, le vecteur X peut aisément être obtenu de l'équation (14) et la solution est unique. C'est en général le cas quand deux points caractéristiques sont observés au moins pendant quatre images consécutives. En utilisant seulement trois images, deux solutions distinctes pour ces angles sont en général fournies, mais dans la plupart des cas, l'une d'entre elles peut être éliminée. Il est également possible de déterminer les angles de roulis et tangage à partir d'un seul point caractéristique, mais dans ce cas, une image supplémentaire est requise. Une fois que le vecteur X est déterminé, la gravité dans le repère absolu RG est aussi déterminée puisque ses composantes sont les trois premières entrées de X. REVENDICATIONS1.- Procédé de détermination, dans un référentiel fixe 3D (RG), de la localisation d'un engin mobile (10) doté d'un capteur d'images (7) associé à un premier référentiel local 3D (R1) correspondant à la prise d'une image, par ledit capteur, de l'environnement de l'engin mobile comprenant au moins deux points caractéristiques dudit environnement, ledit procédé étant caractérisé en ce que les angles de roulis et de tangage de l'engin mobile dans le référentiel fixe correspondant à ladite prise d'image ayant été déterminés dans une étape préalable, les étapes suivantes sont mises en oeuvre, relativement à ladite image prise, pour déterminer l'angle de lacet et/ou les coordonnées de l'engin mobile dans le référentiel fixe : a/ on définit un deuxième référentiel local (Rn) à l'engin mobile comprenant un axe présentant un écart d'angle connu avec l'axe de lacet du référentiel fixe, ledit deuxième référentiel local résultant de la succession de deux rotations, une desdites rotations autour d'un premier axe du premier référentiel local étant fonction de l'angle de roulis (R) déterminé et l'autre desdites rotations étant fonction de l'angle de tangage (P) déterminé ; b/ les coordonnées de chacun des deux points caractéristiques (P1, P2) dans ledit deuxième référentiel local étant égales, à un facteur commun inconnu près, au résultat de la multiplication d'une matrice définissant ladite succession de rotations, par une matrice de paramètres intrinsèques du capteur d'images et par les coordonnées d'un pixel correspondant audit point caractéristique sur l'image, on détermine les valeurs prises par les quotients entre lesdites coordonnées de manière à s'affranchir du facteur commun inconnu ; c/ on exprime ensuite lesdites coordonnées des deux points caractéristiques dans ledit deuxième référentiel local en fonction de la distance (D) entre lesdits points caractéristiques et desdites valeurs de quotients déterminées. 2.- Procédé de détermination de la localisation d'un engin mobile (10) selon la revendication 1, selon lequel on calcule les coordonnées, dans le deuxième référentiel local (Rn), du vecteur joignant les deux points caractéristiques, (P1, P2) et selon lequel l'angle de lacet (YW) de l'engin mobile dans le référentiel fixe correspondant à ladite image prise est déterminé en fonction de la fonction arctan appliquée à une variable fonction desdites coordonnées, dans le repère local, du vecteur joignant les deux points caractéristiques. 3.- Procédé de détermination de la localisation d'un engin mobile (10) selon la revendication 2, selon lequel on calcule les coordonnées de l'engin mobile dans le référentiel fixe (RG) correspondant à ladite image prise, en fonction de coordonnées d'au moins un desdits points caractéristiques et dudit angle de lacet (YW) déterminé dans le référentiel fixe. 4.- Procédé de détermination de la localisation d'un engin mobile selon la revendication 2 ou 3, selon lequel, étant donné trois points caractéristiques (P1, P2, P3) dans le plan perpendiculaire à l'axe de lacet (YW) du référentiel fixe (RG), les étapes suivantes sont mises en oeuvre relativement à ladite image prise : - pour chacune des trois paires de points caractéristiques, on met en oeuvre les étapes b/ et cl et on détermine ainsi trois valeurs d'angles de lacet en fonction des valeurs de roulis et tangage déterminées dans l'étape préalable ; - à partir des trois valeurs d'angles de lacet déterminées, on met en oeuvre une méthode de minimisation d'une fonction cost(1 - (YW 13 YW12 ri)2 +(1w23-dia r2)2 dans laquelle les variables YW13, YW12 et YW23 sont des fonctions non-linéaires de = LRoulis,Tangage]T où yi et 72 représentent deux angles internes du triangle de sommets desdits points caractéristiques ; - en fonction des valeurs de roulis et tangage issues de cette méthode de minimisation, les étapes b/ et cl sont mises en oeuvre relativement à deux des points caractéristiques et on détermine ainsi un angle de lacet dans le référentiel fixe. 5.- Procédé de détermination de la localisation d'un engin mobile (10) selon la revendication 4, comprenant une étape de détermination des angles yi et 72 , selon laquelle on calcule : a c 2 b2 2 a c b 2 2 2 = arccos( ) et 72 - arccos( ) 2ab 2ac où a = \/(/-12 )2 + (y2 - Y1)2 b = ii1113 - 111)2 + (y3 - Y1)2 C = "Ni1113 -112 )2 + (y3 - y2 )2 5 Xi Yi Z. X' les valeurs du/ =- et v =, où z11 sont les coordonnées du point caractéristique Pi dans le deuxième référentiel local, = 1 à 3, étant déterminées à partir 26 ui = Rot.A-1. vi 1 des égalités suivantes : 2. xi Yi , où u,, v, sont les coordonnées d'un pixel Z. de la représentation du point caractéristique P, dans l'image, A est une matrice prédéterminée de paramètres intrinsèques au capteur d'images (7), 2. un facteur de valeur inconnue ; et Rot est la matrice représentant la succession des deux rotations. 6.- Procédé de détermination de la localisation d'un engin mobile (10) selon la revendication 5, selon laquelle ladite étape de détermination des angles yi et 72 est mise en oeuvre dans une phase d'observation d'une méthode de filtre de Kalman. 7.- Procédé de détermination de la localisation d'un engin mobile (10) selon l'une des revendications précédentes, selon lequel le référentiel fixe (RG) est défini tel que le premier point caractéristique en est l'origine, que l'axe de lacet (Z) est de direction contraire au vecteur gravité et que le deuxième point caractéristique est de coordonnée égale à zéro sur l'axe de tangage (Y) du référentiel fixe. 8.- Procédé de détermination de la localisation d'un engin (10) selon l'une des revendications précédentes, selon lequel un axe (Z,-,) du deuxième référentiel local (R,-,) est choisi parallèle à l'axe de lacet (Z) du référentiel fixe (RG). 9.- Procédé de détermination de la localisation d'un engin mobile (10) selon les revendications 7 à 8, selon lequel lorsque les deux points caractéristiques ont une même valeur de coordonnée sur l'axe de lacet (Z) du référentiel fixe, la valeur de coordonnée z de la localisation de l'engin mobile sur l'axe de lacet du référentiel fixe est : z=+ D avec Adu = - ,u2;Av =vi -v2 et D la distance entre les deux 1142 +Av2 points caractéristiques, X2 Y1 ' ret2 = X2 et Y2 - 311 et Y2 sont les coordonnées 111 = V1 = OU V2 = Z1 Z1 Z2 Z2 z1 z2 des deux points caractéristiques dans le deuxième référentiel local (Rn). 10.- Procédé de détermination de la localisation d'un engin mobile (10) selon la revendication 7 et 8, selon lequel lorsque les deux points caractéristiques n'ont pas une même valeur de coordonnée sur l'axe de lacet du référentiel fixe, la valeur de coordonnée z de la localisation de l'engin mobile sur l'axe de lacet (Z) du référentiel fixe est une solution de : (42 ±Av2)( z)2 2(4,12 ±Avv2) Az z ± (A22 avec D D D) 1= 0 D A.,u = pi -#2;4y = y1 -y2 et D la distance entre les deux points caractéristiques, X1 V1 Y1 " X2 ,i. ,, Xi et X2 sont les coordonnées Y2 ,,-, 311 Z1 Y2 Z2 /11 - - I ' 1 - t-v2 = - et v 2 = vu Z1 Z1 Z2 Z2 des deux points caractéristiques dans le deuxième référentiel local. 11.- Procédé de détermination de la localisation d'un engin mobile (10) selon la revendication 9 ou 10, comprenant une étape de détermination de l'angle de lacet YW de l'engin mobile dans le référentiel fixe (RG) à l'aide de la formule suivante : - z(v2 - vi ) YW = -0 = - 2 .arctan , \I(4/12 - /11 ))2 + (Z(/2 - y1 ))2 + (- 4/12 - /11 b 12. Procédé de détermination de la localisation d'un engin mobile (10) selon l'une des revendications 9 à 10, comprenant au moins une des étapes i/ et ii/ ci-dessous, avec - - z(v2 - v1) 0 = - 2 .arctan , .\/(4/12 - /11 »2 ± (Z(/2 - 1/1 )Y + (- 4/12 - /11 b i/ une étape de détermination de la coordonnée x de la localisation de l'engin mobile sur l'axe de roulis du deuxième référentiel fixe (Rn) en calculant x = cos 0. z.,u1 + sin 0. z.v1 ; ii/ une étape de détermination de la coordonnée y de l'engin mobile sur l'axe de tangage (Y) du référentiel fixe en calculant y = -sin 0. z.,u1 + cos 0. z.v1. 13.- Procédé de détermination de la localisation d'un engin mobile (10) selon l'une des revendications 9 à 12, selon lequel /il = 7, y1=3, ,u2= et v2 = A sont 1 z1 z2 Z2 calculées sur la base des égalités suivantes : = Rot.A-1. u1 , où u1, vl sont les coordonnées d'un pixel de la 1 représentation d'un des deux points caractéristiques dans l'image ; et .1. x2 = Rot.A-1 u2 , où u2, v2 sont les coordonnées d'un pixel de la Y2 Z2 V2 1 représentation de l'autre des deux points caractéristiques dans l'image, A est une matrice prédéterminée de paramètres intrinsèques au capteur d'images (7), un facteur de valeur inconnue ; et Rot est la matrice représentant la succession des deux rotations. 14.- Programme d'ordinateur pour localiser, dans un référentiel fixe 3D (RG), un engin mobile (10) doté d'un capteur d'images (7) associé à un premier référentiel local 3D (R1) correspondant à la prise d'une image, par ledit capteur, de l'environnement de l'engin mobile comprenant au moins deux points caractéristiques dudit environnement, ledit programme comportant des instructions pour mettre en oeuvre les étapes d'un procédé selon l'une quelconques des revendications précédentes lors d'une exécution du programme par des moyens de traitement, les angles de roulis et de tangage de l'engin mobile dans le référentiel fixe correspondant à ladite prise d'image ayant été déterminés dans une étape préalable. 15.- Dispositif de détermination, dans un référentiel fixe 3D (RG)' de la localisation d'un engin mobile doté d'un capteur d'images (7) associé à un premier référentiel local 3D (R1) correspondant à la prise d'une image, par ledit capteur, de l'environnement de l'engin mobile comprenant au moins deux points caractéristiques dudit environnement, ledit dispositif de détermination étant adapté pour déterminer préalablement les angles de roulis et de tangage de l'engin mobile dans le référentiel fixe correspondant à ladite prise d'image ; ledit dispositif de détermination étant caractérisé en qu'il est adapté pour déterminer en outre l'angle de lacet et/ou les coordonnées de l'engin mobile dans le référentiel fixe en définissant un deuxième référentiel local (Rn) à l'engin mobile comprenant un axeprésentant un écart d'angle connu avec l'axe de lacet du référentiel fixe, ledit deuxième référentiel local résultant de la succession de deux rotations, une desdites rotations autour d'un premier axe du premier référentiel local étant fonction de l'angle de roulis (R) déterminé et l'autre desdites rotations étant fonction de l'angle de tangage (P) déterminé ; les coordonnées de chacun des deux points caractéristiques (P1, P2) dans ledit deuxième référentiel local étant égales, à un facteur commun inconnu près, au résultat de la multiplication d'une matrice définissant ladite succession de rotations, par une matrice de paramètres intrinsèques du capteur d'images et par les coordonnées d'un pixel correspondant audit point caractéristique sur l'image, ledit dispositif de détermination étant adapté pour déterminer les valeurs prises par les quotients entre lesdites coordonnées de manière à s'affranchir du facteur commun inconnu ; ledit dispositif de détermination étant adapté pour exprimer ensuite lesdites coordonnées des deux points caractéristiques dans ledit deuxième référentiel local en fonction de la distance (D) entre lesdits points caractéristiques et desdites valeurs de quotients déterminées.