FR 3136546 A1 20231215 FR 2205565 A 20220609 Procédé de localisation inertielle mettant en œuvre une corrélation gravimétrique et dispositif associé Un aspect de l’invention concerne un procédé de localisation inertielle d’un porteur utilisant un filtre stochastique de type bayésien récursif mis en œuvre par un dispositif de localisation inertielle, le procédé comprenant une pluralité de cycles, chaque cycle comportant une étape de propagation et, lorsqu’un mesure est disponible, une étape de mise à jour, procédé dans lequel, lors de l’étape de propagation, la propagation est effectuée à l’aide d’un modèle de propagation prenant en compte au moins une composante d’une matrice des gradients spatiaux du vecteur de pesanteur fournie à l’aide par un modèle gravimétrique enrichi. Figure 1 Procédé de localisation inertielle mettant en œuvre une corrélation gravimétrique et dispositif associé DOMAINE TECHNIQUE DE L’INVENTION Le domaine technique de l’invention est celui de la localisation inertielle. La présente invention concerne un procédé de localisation inertielle et en particulier un procédé de localisation inertielle par corrélation gravimétrique par propagation d’erreur. ARRIERE-PLAN TECHNOLOGIQUE DE L’INVENTION L’estimation de la position d’un mobile sans aide GNSS peut être réalisée par un équipement constitué d’un gravimètre et d’une centrale de localisation inertielle. Une série de couples (position, gravité) est recherchée dans une carte gravimétrique représentative de la zone traversée d’après les mesures du gravimètre et la position estimée par la centrale pour corriger la position, où le signal gravimétrique est une fonctionnelle gravimétrique telle que la pesanteur, la gravitation, l’anomalie de gravité, la perturbation de gravité. Plusieurs méthodes existent pour faire la recherche des couples (position, signal gravimétrique) dans la carte. Une première méthode repose sur un modèle simplifié de l’erreur de position. Sur un intervalle de temps donné, cette erreur est modélisée comme une constante inconnue. Une estimation de cette constante est faite en décalant la trajectoire estimée dans la zone de recherche et en recherchant un pic de corrélation entre les valeurs lues sur la carte le long de la trajectoire décalée et les mesures relevées sur cet intervalle de temps. Cette méthode est particulièrement utile lorsque l’erreur de position est de l’ordre de 10 km ou plus. Une deuxième méthode utilisant la carte gravimétrique repose sur une hybridation serrée entre la centrale de localisation inertielle et le gravimètre. L’état de navigation est alors estimé de façon optimale à chaque observation d’après toutes les mesures inertielles et observations précédentes. Ces deux méthodes ont un problème en commun : le fait que le gravimètre soit en mouvement dégrade la précision de la mesure et limite donc la correction de position. Plus la vitesse est élevée, moins la position est corrigée. Cela s’explique en partie par le temps de mesure du gravimètre et les variations du champ auquel est soumis le gravimètre durant la mesure. Un temps de mesure de 10 secondes à la vitesse de 900 km/h peut conduire à des variations de 10 µg pendant la mesure dans certaines zones géographiques alors que la précision intrinsèque d’un gravimètre absolu est typiquement de l’ordre de 10 ng. Le coût de l’équipement devient donc élevé par rapport à la qualité de la mesure quand la vitesse augmente. Le GNSS et le gravimètre sont des moyens de mesure passifs présentant un avantage de discrétion, mais le GNSS n’est pas toujours disponible et le gravimètre est moins performant à grande vitesse. Aussi, il existe un besoin d’un procédé de localisation inertielle sans contrainte de disponibilité qui peut être mis en œuvre sur des porteurs se déplaçant à grande vitesse. L’invention offre une solution aux problèmes évoqués précédemment, en proposant un procédé dans lequel le gravimètre extérieur peut être remplacé par un capteur dont la mesure est une fonction de la position ou de la vitesse du porteur, et dans lequel un modèle gravimétrique enrichi est appliqué dans le modèle de propagation de l’état de navigation (position, vitesse, attitude) à chaque position de la zone d’incertitude modélisée par la centrale, le terme « enrichi » signifiant que ce modèle permet de décrire le champ de pesanteur d’une façon plus fine que le modèle normal utilisé par défaut dans les centrales de localisation inertielle, le modèle normal de pesanteur étant un terme utilisé en géodésie physique désignant le champ de pesanteur généré par un ellipsoïde de référence équipotentiel en rotation modélisant la surface d’une planète et qui, dans l’état de l’art concernant la planète Terre, est celui décrit par le système de référence WGS84 ou GRS80. Pour cela, un premier aspect de l’invention concerne un procédé de localisation inertielle d’un porteur utilisant un filtre stochastique de type bayésien récursif mis en œuvre par un dispositif de localisation inertielle, le procédé comprenant une pluralité de cycles, chaque cycle comportant une étape de propagation et une étape de mise à jour, procédé dans lequel, lors de l’étape de propagation, la propagation est effectuée à l’aide d’un modèle de propagation prenant en compte au moins une composante d’une matrice des gradients spatiaux du vecteur de pesanteur fournie par un modèle gravimétrique enrichi. Ainsi, lors de la propagation, l’état de navigation (position, vitesse, attitude) est calculé d’après l’état de navigation initial de chaque cycle, d’après les mesures inertielles et d’après le modèle gravimétrique enrichi. Le procédé selon l’invention est particulièrement avantageux en ce que la lecture du modèle gravimétrique enrichi à la position estimée, différente de la position vraie, crée, en raison de la répartition spatiale non-uniforme du vecteur gravimétrique enrichi, une corrélation entre l’erreur de vitesse et l’erreur de position dans le voisinage de la position estimée, et en ce que le filtre stochastique calcule l’incertitude en modélisant cette corrélation lors de l’étape de propagation par linéarisation du vecteur gravimétrique enrichi fonction de la position à chaque instant de propagation sous forme d’une matrice de gradients spatiaux utilisée partiellement. Ainsi, une observation partielle de la vitesse permet alors de corriger une partie de la position par l’intermédiaire de cette corrélation lors de l’étape de mise à jour. De façon similaire, une observation partielle de la position permet de corriger une partie de la vitesse et une autre partie de la position, lors de cette même étape. Grâce à l’invention, il est donc possible de mettre en œuvre un procédé de localisation sur des porteurs se déplaçant à grande vitesse en modélisant la propagation de l’erreur sur des intervalles de temps inférieurs au cycle associé au filtre stochastique (par exemple au cycle de Kalman lorsqu’un tel filtre est utilisé). De plus, ce procédé peut être mis en œuvre à l’aide d’une mesure non susceptible d’être détectée, par exemple une mesure de l’altitude du porteur ou bien encore de sa vitesse verticale. En outre, sa mise en œuvre ne sollicite que peu de ressources de calcul. Outre les caractéristiques qui viennent d’être évoquées dans le paragraphe précédent, le procédé selon un premier aspect de l’invention peut présenter une ou plusieurs caractéristiques complémentaires parmi les suivantes, considérées individuellement ou selon toutes les combinaisons techniquement possibles. Dans un mode de réalisation le filtre bayésien récursif linéarise la loi de propagation. Dans un mode de réalisation, le dispositif de localisation inertielle comprend un moyen de calcul, un bloc de capteurs inertiels et au moins un capteur externe, et pour chaque cycle avec un entier positif non nul : lors de l’étape de propagation, l’état de navigation en fin de cycle est calculé en fonction de l’état de navigation en début de cycle , de mesures inertielles, et du vecteur de pesanteur obtenu d’après le modèle gravimétrique enrichi, l’état d’erreur est calculé en fonction de l’état d’erreur , et la matrice de covariance est calculée en fonction de la matrice de covariance , le calcul de l’état d’erreur et de sa covariance s’effectuant à l’aide d’une matrice de transition au cycle et d’une matrice de covariance de bruit de modèle au cycle de sorte que : la matrice de transition étant déterminée à l’aide d’au moins une composante de la matrice des gradients spatiaux du vecteur de pesanteur obtenus d’après le modèle gravimétrique enrichi, une matrice d’observation étant ensuite déterminée à partir de l’état de navigation en fin de cycle et d’un modèle d’observation ; lors de l’étape de mise à jour, l’état de l’erreur après mis à jour et la matrice de covariance de l’erreur d’estimation après mise à jour sont calculés, par le moyen de calcul et à l’aide du filtre stochastique, à partir de l’état de l’erreur avant mise à jour , de la matrice de covariance de l’erreur d’estimation avant mise à jour , de la matrice d’observation , d’une matrice de covariance du bruit de mesure , et d’une observation relative à au moins une fonction de la vitesse et/ou de la position du porteur effectuée par le capteur externe au cours ou à la fin de l’étape de propagation de sorte à réduire au moins en partie l’état d’erreur après mise à jour. Dans un mode de réalisation, chaque cycle de filtrage k est divisé en intervalles de temps ainsi qu’en P intervalles de temps de sorte que avec où et et sont des entiers positifs non nuls, où est la période du filtre, l’étape de propagation de chaque cycle de filtrage comprenant, à partir d’une matrice de propagation : Pour chaque intervalle repéré par l’indice avec compris entre 1 et N, une sous-étape de calcul de l’état de navigation par propagation de l’état à partir des mesures du bloc de capteurs inertiels, du modèle gravimétrique enrichi à la position de l’état ; Pour chaque intervalle repéré par l’indice compris entre , et , étant un entier entre 1 et P désignant l’indice d’un intervalle avec où est la partie entière de x, une sous-étape de calcul d’une matrice de transition élémentaire initialisée à la matrice identité au début de l’intervalle de temps repéré par l’indice et complétée par intégration de la matrice de propagation par rapport au temps sur les intervalles de temps composant l’intervalle , la matrice obtenue par l’intégration de la matrice de propagation sur un intervalle de temps étant ajoutée à la matrice obtenue lors de l’intégration sur l’intervalle de temps précédent de sorte à constituer progressivement cette matrice de transition élémentaire , et une matrice de transition sur le cycle k initialisée à la matrice identité étant calculée progressivement à chaque nouvelle valeur de par le produit matriciel ; A l’issue du dernier intervalle du cycle k et donc du dernier intervalle ( et ), une sous-étape de détermination de la matrice d’observation à partir de l’état de navigation propagé et d’un modèle d’observation, et une sous-étape de propagation de l’état d’erreur et de la matrice de covariance de cet état : avec la matrice de transition obtenue à la fin du dernier intervalle au cycle considéré. Dans un mode de réalisation, le dispositif de localisation inertielle comprend un télémètre laser et/ou un lidar, et une mesure de la position et/ou de la vitesse du porteur est réalisée par télémétrie lors de l’étape de mise à jour. Dans un mode de réalisation, le dispositif de localisation comprend une caméra, et une mesure de la position et/ou de la vitesse du porteur est réalisée par caméra lors de l’étape de mise à jour. Dans un mode de réalisation, le dispositif de localisation comprend un moyen de mesure de la position verticale, et une mesure de la position et/ou de la vitesse verticale du porteur est réalisée par ledit moyen de mesure lors de l’étape de mise à jour. Dans un mode de réalisation, le dispositif de localisation inertielle est à composants liés (dit strapdown en anglais) mécanisée dans le repère géographique local (Nord, Ouest, Haut) ou (Nord, Est, Bas). Dans un mode de réalisation, le dispositif de localisation inertielle est à composants liés mécanisée dans un repère dont un des axes coïncide avec l’axe vertical du repère géographique local, dit repère en azimut libre. Ce repère permet notamment de passer les pôles géographiques sans provoquer de division par zéro dans les calculs. Dans un mode de réalisation, le dispositif de localisation inertielle est à composants liés mécanisée dans le repère géocentrique lié à la planète (en général noté [t]). Dans un mode de réalisation, le dispositif de localisation inertielle est à composants liés mécanisée dans le repère géocentrique inertiel (en général noté [i]). Dans un mode de réalisation, le dispositif de localisation inertielle est un dispositif à cardans. Un deuxième aspect de l’invention concerne un dispositif de localisation inertielle comprenant un bloc de capteurs inertiels, au moins un capteur externe, une mémoire de stockage du modèle gravimétrique enrichi, et des moyens configurés pour mettre en œuvre un procédé selon un premier aspect de l’invention. Un troisième aspect de l’invention concerne un programme d’ordinateur comprenant des instructions qui conduisent le dispositif selon un deuxième aspect de l’invention, lorsque ces instructions sont exécutées par le dispositif, à mettre en œuvre le procédé selon un premier aspect de l’invention. Un quatrième aspect de l’invention concerne un support lisible par ordinateur, sur lequel est enregistré le programme d’ordinateur selon un troisième aspect de l’invention. L’invention et ses différentes applications seront mieux comprises à la lecture de la description qui suit et à l’examen des figures qui l’accompagnent. BREVE DESCRIPTION DES FIGURES Les figures sont présentées à titre indicatif et nullement limitatif de l’invention. La montre une représentation schématique du principe de corrélation gravimétrique selon l’invention. La montre une représentation schématique de l’utilisation d’un filtre de Kalman pour la mise en œuvre de la corrélation gravimétrique. La montre un logigramme d’un procédé de localisation selon l’invention. La illustre le fonctionnement de l’algorithme de calcul de la solution (position, vitesse, attitude) en fonction des mesures inertielles. La montre une représentation schématique d’un dispositif de localisation selon l’invention. La compare le comportement de l’erreur sur la position horizontale et le cap, pour des capteurs inertiels isolés de l’environnement vibratoire et thermique, entre l’état de la technique et le procédé selon l’invention. La illustre le chronogramme des calculs lors de l’étape de propagation. La illustre la différence entre états en boucle fermée et états en boucle ouverte. La illustre la définition du repère géocentrique [t] lié à la planète et du repère géographique local (Nord, Ouest, Haut) [g] ainsi que la définition de la latitude L, longitude G et altitude ellipsoïdale Z appelée aussi . La illustre la définition du repère géocentrique [t] lié à la planète et du repère géocentrique inertiel [i]. DESCRIPTION DETAILLEE Sauf précision contraire, un même élément apparaissant sur des figures différentes présente une référence unique. Introduction Les méthodes de localisation inertielle sont en général mises en œuvre par une centrale de localisation inertielle qui regroupe un bloc de capteurs inertiels (ou BSI comprenant des accéléromètres et des gyroscopes ou des gyromètres), un moyen de calcul et, de préférence, des éléments d’amortissement mécaniques sur lesquels est habituellement monté le bloc de capteurs inertiels afin de minimiser l’impact des chocs et vibrations sur la précision de localisation. Cette centrale de localisation inertielle peut fonctionner dans un premier mode, dit mode d’alignement, qui correspond à l’initialisation de l’état de navigation à l’aide des mesures du bloc de capteurs inertiels et/ou d’un capteur externe. La présente invention ne concerne pas ce premier mode. Quand l’alignement est terminé, la centrale de localisation peut passer en mode navigation : c’est ce mode de navigation qui est l’objet de l’invention. Dans ce mode, les estimations de l’état de navigation sont calculées selon des cycles à l’aide d’un filtre stochastique (appelé aussi filtre de navigation), chaque cycle comportant aux moins deux étapes : une étape de propagation de l’état de navigation et de l’état du filtre stochastique ; et une étape de mise à jour de l’état de navigation et de l’état du filtre stochastique. Le filtre stochastique le plus utilisé en navigation est le filtre de Kalman étendu. Le calcul de propagation est partagé entre un algorithme de propagation non-linéaire appelé « navigateur » et un algorithme de propagation d’erreurs avec des lois linéarisées d’après les estimations du navigateur à différents instants. Ce dernier algorithme considère que les erreurs sont gaussiennes et calcule l’espérance mathématique et la matrice de covariance des erreurs propagées. La prise en compte des observations est effectuée par le filtre de Kalman étendu dans une étape de mise à jour. C’est le fonctionnement standard d’une centrale de localisation inertielle. Le procédé selon l’invention se distingue du fonctionnement standard décrit ci-dessus en ce que le modèle de propagation utilisé lors de l’étape de propagation pour propager l’état d’erreur noté dans la suite prend en compte au moins une composante, de préférence une pluralité de composantes, de la matrice des gradients spatiaux du vecteur gravimétrique d’un modèle gravimétrique enrichi, ce ou ces composantes étant lues dans le modèle à la position du porteur estimée par le navigateur. Dans un mode de réalisation, on entend par modèle gravimétrique enrichi un modèle dont la résolution du modèle est meilleure que 15 minutes d’arc et l’erreur RMS est inférieure à 15 mGal pour la plupart des zones géographiques (80% en surface ou plus) susceptibles d’être parcourues par le porteur. La résolution est en général définie comme le pas de grille moyen des mesures ramenées à l’ellipsoïde de référence ayant permis de construire le modèle enrichi. Pour un modèle global décomposé en harmoniques sphériques, cela correspond à un modèle dont le plus grand degré d’harmonique sphérique est supérieur à 700 environ. Dans un mode de réalisation préférentiel, les erreurs accélérométriques fournies par le bloc de capteurs inertiels sont inférieures à quelques 10 mGal environ. Les filtres stochastiques Bien que, dans la suite, l’invention soit illustrée à l’aide d’un filtre de Kalman étendu, elle ne se cantonne pas à ce seul filtre stochastique. Toutes les variantes du filtre de Kalman étendu dans lesquelles l’étape de propagation est réalisée avec une loi de propagation linéarisée peuvent être exploitées : par exemple des filtres améliorant la vitesse de convergence, ou des filtres améliorant des problèmes numériques, ou des filtres améliorant des problèmes de robustesse. Bien entendu, il existe encore d’autres variantes de filtres stochastique, mais le principe de l’invention (c’est-à-dire la corrélation gravimétrique par propagation) reste toujours valide bien qu’il puisse prendre des formes différentes. Modèles gravimétriques enrichis Comme mentionné précédemment, dans le procédé selon l’invention, le modèle gravimétrique est un modèle gravimétrique enrichi permettant de définir le gradient spatial du vecteur de pesanteur à la position estimée du porteur. Par contraste, dans les centrales de localisation inertielle de l’état de la technique, le modèle de pesanteur utilisé est le modèle normal ou une approximation dudit modèle, ce dernier comportant 4 paramètres géodésiques qui, dans le cas de la planète Terre sont définis dans des systèmes de référence tels que le WGS84 ou le GRS80. Le modèle normal décrit le champ de pesanteur généré par une surface lisse ellipsoïdale, équipotentielle, en rotation, dotée de la masse de la planète. L’effet du relief topographique sur le champ de pesanteur n’est pas modélisé. L’écart entre le champ de pesanteur vrai et le champ de pesanteur normale à une même position est appelé champ perturbateur de gravité. Il se traduit par des déviations de verticale et des erreurs sur le module de pesanteur. Les déviations de verticale excitent les erreurs de navigation et produisent des oscillations de Schuler en position et vitesse horizontales pouvant potentiellement diverger mais pouvant être partiellement amorties par des observations. Certes, dans l’état de la technique, des modèles de déviation de verticale sont employés pour réduire ces oscillations dans différents systèmes de navigation. Ces modèles sont par définition plus riches que le modèle normal. D’autre part, des gradients spatiaux du modèle normal sont modélisés dans le modèle de propagation des erreurs du filtre de Kalman étendu. Cependant, dans l’état de la technique, les gradients spatiaux du vecteur de pesanteur enrichi ne sont pas modélisés dans la matrice de propagation. Il en résulte que, dans les systèmes de navigation utilisant un modèle gravimétrique enrichi, une corrélation non modélisée entre l’erreur de vitesse et l’erreur de position apparaît due au fait que, dans le navigateur, le modèle gravimétrique enrichi est lu à la position estimée et non à la position vraie. Outre cette incohérence, il apparaît que l’erreur de vitesse contient une information sur l’erreur de position, et que certains capteurs existants peuvent permettre de capter cette information afin d’estimer l’erreur de positon. La performance atteignable dépend de la richesse du modèle gravimétrique et de la qualité des accéléromètres. Des modèles mondiaux parmi lesquels l’EGM2008, l’EIGEN-6C4, l’XGM2019, construits d’après des mesures terrestres et spatiales du champ de pesanteur permettent de modéliser le champ de pesanteur avec une résolution de l’ordre de 10 km. D’autres modèles mondiaux de l’état de l’art construits d’après des modèles de topographie réduisent cette résolution à 6 km. De meilleures résolutions sont obtenues par certains modèles locaux. Ces modèles comportent des erreurs modélisées partiellement. Illustration schématique du p rincipe de l’invention Avant d’illustrer en détail l’invention, il peut être utile de présenter de manière simplifiée, dans un mode de réalisation où le filtre utilisé est un filtre de Kalman étendu, le principe sur lequel elle repose (c’est-à-dire la corrélation gravimétrique par propagation d’erreur) afin de voir comment, à partir d’au moins une composante du gradient spatial du vecteur de pesanteur enrichi, il est possible d’améliorer l’évaluation de la position et/ou de la vitesse et/ou de l’attitude en réduisant l’erreur sur ces dernières. Cette présentation permettra également de montrer comment le procédé selon l’invention peut être mis en œuvre même en l’absence de gravimètre, un capteur de position ou de vitesse étant suffisant. Dans l’exemple simplifié décrit dans la suite, il est supposé que le modèle de pesanteur est parfait ou bien encore que les mesures des accéléromètres sont parfaites. Bien entendu, il n’en est rien dans la réalité et les opérations présentées dans les paragraphes qui suivent doivent être faites en tenant compte des erreurs sur le modèle et/ou les mesures. Sur la sont représentées deux positions d’un porteur, la position réelle PR du porteur, c’est-à-dire la position où le porteur se situe réellement, et la position estimée PE de ce porteur, c’est-à-dire la position du porteur estimée par une centrale de localisation inertielle présente dans le porteur. Les accéléromètres de la centrale de localisation inertielle mesurent une force spécifique FS à la position réelle PR. Pour obtenir l’accélération du porteur, le navigateur complète cette mesure par le vecteur de pesanteur ou de gravitation (le vecteur de pesanteur est la somme du vecteur de gravitation correspondant à l’attraction universelle et à l’accélération centrifuge due à la rotation de la planète sur laquelle se trouve le porteur), ainsi que des accélérations dépendant du type de mécanisation utilisé par le navigateur et connues de l’état de l’art (la mécanisation définit le repère dans lequel est exprimée la navigation). Ce vecteur de pesanteur est issu d’un modèle de pesanteur enrichi et est calculé à la position estimée PE. Cette accélération permet alors de calculer la vitesse du porteur à l’instant suivant et donc également sa position. Aussi, bien que la force spécifique FS soit issue d’une mesure, la valeur de l’accélération du porteur est fonction d’une valeur estimée du vecteur de pesanteur calculée à l’aide d’un modèle enrichi à partir de la position estimée PE. Aussi, une erreur sur la position estimée PE se propage, par l’intermédiaire de la prise en compte de la pesanteur, en une erreur sur l’accélération du porteur (et donc sur la vitesse et la position dudit porteur, ces grandeurs étant obtenues par intégration). Sur la , l’erreur de position sur l’axe est notée . Cette erreur de position sur l’axe engendre une différence entre la pesanteur selon l’axe à la position réelle PR et la pesanteur déterminée à l’aide d’un modèle à la position estimé PE, noté sur la . Ces deux erreurs peuvent être liées l’une à l’autre à l’aide de la relation suivante : Où est l’erreur sur la composante en de la pesanteur, est le gradient spatial selon l’axe de la composante en de la pesanteur à la position estimée PE et est l’erreur sur la composante en de la position. Comme déjà mentionné, cette erreur sur la valeur de la pesanteur selon l’axe se traduit directement par une erreur dans l’accélération selon ce même axe selon la relation : Il est possible de dériver des relations similaires pour l’erreur sur la composante en de la vitesse de sorte que, pendant un intervalle de temps T pendant lequel le gradient peut être considéré constant, la contribution de l’erreur de position à l’erreur de vitesse z, hormis tous les autres contributeurs, est : Autrement dit, il est possible de propager l’erreur sur l’erreur en établissant une corrélation (gravimétrique) entre ces deux grandeurs. En outre, cette propagation ne nécessite pas de mesure externe du vecteur de pesanteur ou de son gradient spatial (en général, des états du filtre stochastique sont réservés à l’estimation d’une correction du modèle de pesanteur, cette estimation pouvant être considérée comme une mesure interne du vecteur de pesanteur), la valeur du vecteur et du gradient spatial pouvant par exemple être fournies par un modèle gravimétrique enrichi. Une observation sur la voie verticale permet alors de corriger la position horizontale. La relation précédente constitue évidemment une approximation servant à illustrer l’invention et d’autres termes additionnels sont présents de sorte que la relation de proportionnalité [Math. 3] n’est pas vérifiée. Cependant, cette approximation permet de présenter de manière simplifiée le principe de l’invention. Procédé de localisation inertielle Comme illustré à la , à la , à la et à la , un premier aspect de l’invention concerne un procédé 100 de localisation inertielle d’un porteur à l’aide d’un filtre stochastique mis en œuvre par un dispositif DI de localisation inertielle. Le procédé 100 selon l’invention comprend une pluralité de cycles, chaque cycle comprenant une étape 1E1 de propagation et une étape 1E2 de mise à jour. De plus, dans le procédé 100 selon l’invention, lors de l’étape 1E1 de propagation, la propagation est effectuée à l’aide d’un modèle de propagation prenant en compte au moins une composante de la matrice des gradients spatiaux du vecteur de pesanteur fournie à l’aide d’un modèle gravimétrique enrichi. Dans un mode de réalisation, le filtre stochastique est un filtre de Kalman étendu. Bien entendu, comme détaillé en introduction, d’autres filtres peuvent être utilisés dans le cadre de l’invention, mais le filtre de Kalman étendu a été choisi pour illustrer l’invention car il s’agit du filtre le plus courant dans le domaine de la localisation. Dans un mode de réalisation illustré à la , le dispositif DI de localisation inertielle comprend un moyen de calcul MC (par exemple un processeur ou une carte ASIC), un bloc de capteurs inertiels BSI (ou ISB de l’anglais Inertial Sensor Block), au moins un capteur externe CE, et une mémoire de stockage MS configurée pour stocker les instructions et les données nécessaires à la mise en œuvre d’un procédé 100 selon l’invention, et notamment le modèle gravimétrique enrichi. Dans un mode de réalisation préférentiel, un moyen de stockage spécifique est utilisé pour stocker le modèle gravimétrique enrichi, par exemple une deuxième mémoire de stockage distincte de celle stockant les instructions de mise en œuvre du procédé 100 selon l’invention. La sera utilisée pour illustrer les différentes étapes 1E1, 1E2 du procédé 100 selon l’invention. Dans cette , l’invention est illustrée sur deux axes (l’axe et l’axe ), et une mesure (ou observation) de la vitesse selon l’axe . Ainsi, dans cette illustration simplifiée, l’état d’erreur qui, dans en réalité, contient non seulement les erreurs de navigation mais aussi des erreurs d’états complémentaires contribuant à réduire les erreurs de navigation, est réduit à l’erreur sur le vecteur d’état notée à l’instant et décrite par : Où est l’incertitude sur la composante en , est l’incertitude sur la composante en et est l’incertitude sur la composante en de la vitesse. De plus, toujours afin d’illustrer le principe de l’invention dans le cas d’une modélisation gaussienne des erreurs, sur la , les ellipsoïdes d’incertitude à 1σ sont représentées en projection dans trois plans : Le plan ( en bas à gauche ; Le plan (Xu 1 , ) en haut à gauche où Xu 1 représente l’axe du gradient spatial du vecteur de pesanteur de direction ; Le plan ( , Xu 2 ) en bas à droite où Xu 2 représente l’axe orthogonal au gradient. Comme évoqué en introduction, le procédé 100 selon l’invention comprend deux étapes : une étape 1E1 de propagation et une étape 1E2 de mise à jour, ces deux étapes se répétant selon une période , chaque répétition correspondant à un cycle avec un entier naturel non nul. La période est la période associée au filtre stochastique, ici le filtre de Kalman étendu. Ces deux étapes vont être détaillées dans la suite. Dans un mode de réalisation, le modèle gravimétrique enrichi comprend un ensemble de données et de fonctions (par exemple des interpolations ou des fonctions d’harmoniques sphériques ou ellipsoïdales) fournissant un vecteur gravimétrique et une fonction de la matrice de gradients spatiaux de ce vecteur fournissant au moins un gradient spatial permettant de réaliser des approximations telles que des valeurs de gradients moyennes locales, en fonction de la position du porteur. Par exemple, dans le cas d’une grille, le taux de variation spatiale entre deux points de grille est aussi une approximation du gradient dans cette direction : . Ce taux de variation est donc la valeur moyenne du gradient spatial de dans la direction . Dans un mode de réalisation alternatif, le modèle gravimétrique enrichi fournit un vecteur gravimétrique et une matrice de gradients spatiaux de ce vecteur correspondant au vecteur de pesanteur et une fonction de son jacobien dans un repère [g] (ou repère géographique) à des positions exprimées en coordonnées géographiques relativement à l’ellipsoïde de référence où L est la latitude géographique, G la longitude et l’altitude ellipsoïdale, où est une fonction d’approximation du jacobien permettant, par exemple, de forcer à zéro certains termes ou réaliser une moyenne locale d’après des paramètres, ainsi qu’un modèle d’erreur associé, le vecteur gravimétrique et sa matrice de gradients spatiaux étant et où [q] représente le repère de mécanisation correspondant au repère sur lequel est asservie la localisation inertielle. Par exemple le repère [g] rendant le repère plateforme de la localisation inertielle horizontal tourné vers le Nord. Si la localisation inertielle est à cardans, le BSI est fixe par rapport à une plateforme mécanique rendue horizontale tournée vers le Nord. Si la localisation est de type strapdown, le BSI est fixe par rapport au porteur et une plateforme virtuelle est calculée, équivalente à la plateforme mécanique de la centrale à cardans. Dans l’état de l’art, les centrales sont toutes de type strapdown. Dans certains cas, la centrale strapdown peut elle-même être montée sur une plateforme à cardans, rendant la plateforme virtuelle fixe par rapport à la plateforme physique. Pour mémoire, le repère plateforme en azimut libre connu de l’homme du métier permet de passer les pôles géographiques et correspond à un repère déduit de [g] par un angle autour de l’axe où est entretenu par la vitesse Ouest du porteur, ou par exemple le repère terrestre géocentrique [t]. De plus : [Math. 5] Où sont les coordonnées du vecteur de pesanteur dans le repère [g] à la position définie par la latitude L, la longitude G et l’altitude ellipsoïdale , et où une composante telle que représente la variation de lorsque la position linéaire (en mètres par opposition à la position angulaire comme L et G) évolue le long de l’axe y de [g] de la valeur ., le vecteur de pesanteur enrichi et sa matrice de gradients spatiaux étant calculés par : Dans un mode de réalisation alternatif, le modèle gravimétrique enrichi fournit un vecteur gravimétrique et une matrice de gradients spatiaux de ce vecteur correspondant au vecteur perturbateur de gravité et une fonction de son jacobien dans le repère [g] à des positions exprimées en coordonnées géographiques relativement à l’ellipsoïde de référence où L est la latitude géographique, G la longitude et l’altitude ellipsoïdale, où est une fonction d’approximation du jacobien permettant par exemple de forcer à zéro certains termes ou réaliser une moyenne locale d’après des paramètres, ainsi qu’un modèle d’erreur associé, le vecteur gravimétrique et sa matrice de gradients spatiaux étant et avec : Où sont les coordonnées du vecteur perturbateur de gravité dans le repère [g] à la position définie par la latitude L, la longitude G et l’altitude ellipsoïdale , et où une composante telle que représente la variation de lorsque la position linéaire (par opposition à la position angulaire) évolue le long de l’axe y de [g] de la valeur . , le vecteur de pesanteur enrichi et sa matrice de gradients spatiaux étant calculés par : où , sont définis à la même altitude et à la même position horizontale, le calcul du vecteur de pesanteur enrichi pouvant être réalisé de façon équivalente dans l’algorithme de calcul de l’état de navigation en intégrant la vitesse avec le vecteur et en complétant les mesures accélérométriques à l’aide de , une approximation possible, valable jusqu’à 100 km d’altitude pour la planète Terre, étant dans le repère (Nord, Ouest, Haut – un tel repère est illustré à la ) : Avec : Où sont les paramètres géodésiques du modèle de pesanteur normale représentant respectivement le rayon équatorial, l’aplatissement, la vitesse angulaire de rotation dans l’espace inertiel, et le produit de la constante gravitationnelle universelle par la masse de la planète. Dans un mode de réalisation alternatif, le modèle gravimétrique enrichi fournit un vecteur gravimétrique et une matrice de gradients spatiaux de ce vecteur correspondant au vecteur d’anomalie de gravité et une fonction de son jacobien dans le repère [g] à des positions exprimées en coordonnées géographiques relativement à l’ellipsoïde de référence où L est la latitude géographique, G la longitude et l’altitude ellipsoïdale, où est une fonction d’approximation du jacobien permettant par exemple de forcer à zéro certains termes ou réaliser une moyenne locale d’après des paramètres, ainsi qu’un modèle d’erreur associé, le vecteur gravimétrique et sa matrice de gradients spatiaux étant et avec : Où sont les coordonnées du vecteur d’anomalie de gravité dans le repère [g] à la position définie par la latitude L, la longitude G et l’altitude ellipsoïdale , et où une composante telle que représente la variation de lorsque la position linéaire (par opposition à la position angulaire) évolue le long de l’axe y de [g] de la valeur , le vecteur de pesanteur enrichi et sa matrice de gradients spatiaux étant calculés par : Où sont définis à l’altitude ellipsoïdale à la position horizontale (L,G), et où sont définis à la surface de l’eau à la même position horizontale (L,G). Ce cas de figure est réservé à la navigation maritime. Dans ce cas, peuvent être calculés de façon plus simple que [Math. 9] et [Math. 10] avec les expressions suivantes : le calcul du vecteur de pesanteur enrichi pouvant être réalisé de façon équivalente dans l’algorithme de calcul de l’état de navigation en intégrant la vitesse avec le vecteur et en complétant les mesures accélérométriques à l’aide de . Dans un mode de réalisation, le modèle gravimétrique enrichi fournit un vecteur gravimétrique et une matrice de gradients spatiaux de ce vecteur correspondant au vecteur de gravitation , ce dernier étant l’attraction exercée par la masse de la planète, et une fonction de son jacobien dans le repère [g] à des positions exprimées en coordonnées géographiques relativement à l’ellipsoïde de référence où L est la latitude géographique, G la longitude et l’altitude ellipsoïdale, où est une fonction d’approximation du jacobien permettant par exemple de forcer à zéro certains termes ou réaliser une moyenne locale d’après des paramètres, ainsi qu’un modèle d’erreur associé, le vecteur gravimétrique et sa matrice de gradients spatiaux étant et avec : Où sont les coordonnées du vecteur de gravitation dans le repère [g] à la position définie par la latitude L, la longitude G et l’altitude ellipsoïdale , et où une composante telle que représente la variation de lorsque la position linéaire (par opposition à la position angulaire) évolue le long de l’axe y de [g] de la valeur , le vecteur de pesanteur enrichi et sa matrice de gradients spatiaux étant calculés par : Où est l’accélération centrifuge due à la rotation de la planète dans l’espace et où sont calculés à la même position (L,G, inertiel, avec : Bien entendu, certains termes peuvent être approximés lors du calcul. Dans un mode de réalisation alternatif, le modèle gravimétrique enrichi fournit un vecteur gravimétrique et une matrice de gradients spatiaux de ce vecteur correspondant au vecteur de gravitation et une fonction de son jacobien dans le repère géocentrique [t] (un tel repère est illustré à la lié à la planète à des positions exprimées en coordonnées cartésiennes ), où est une fonction d’approximation du jacobien permettant par exemple de forcer à zéro certains termes ou réaliser une moyenne locale d’après des paramètres, ainsi qu’un modèle d’erreur associé, le vecteur gravimétrique et sa matrice de gradients spatiaux étant et avec : Où sont les coordonnées du vecteur de gravitation dans le repère [t] à la position définie par la latitude L, la longitude G et l’altitude ellipsoïdale , et où une composante telle que représente la variation de lorsque la position linéaire (par opposition à la position angulaire) évolue le long de l’axe y de [t] de la valeur , le vecteur de pesanteur enrichi et sa matrice de gradients spatiaux étant calculés par : Où sont calculés à la même position (L,G, correspondant aux coordonnées cartésiennes ( dans [t] avec : D’autres modèles enrichis peuvent également être utilisés dans le cadre de l’invention. L’é tape 1E1 d e propagation Le procédé 100 selon l’invention comprend une étape 1E1 de propagation lors de laquelle l’état de navigation en fin de cycle est calculé en fonction de l’état de navigation de début de cycle , de mesures inertielles, du vecteur de pesanteur obtenu d’après le modèle de pesanteur enrichi (et des lois de la dynamique de Newton : cette propagation est effectuée selon des équations connues du domaine et détaillées en partie dans la suite, les équations en question dépendant du repère de mécanisation choisi). De même, lors de cette étape, l’état d’erreur en fin de propagation , est calculé en fonction de l’état d’erreur en début de cycle , d’une matrice de transition , et la matrice de covariance en fin de propagation est calculée en fonction de la matrice de covariance en début de cycle , de la matrice de transition , d’une matrice de covariance de bruit de modèle , de sorte que : De plus, la matrice de transition est déterminée à l’aide d’au moins une composante, de préférence une pluralité de composantes, de la matrice des gradients spatiaux du vecteur de pesanteur donné par le modèle de pesanteur enrichi. Autrement dit, lors de cette étape 1E1, au cycle , l’état d’erreur et sa matrice de covariance en début de cycle sont propagés en aboutissant à l’état d’erreur et sa matrice de covariance en fin de cycle . De plus, lors de cette étape 1E1, une matrice d’observation est déterminée à partir de l’état de navigation propagé et d’un modèle d’observation. Cette étape 1E1 de propagation est illustrée à la où l’erreur en début de cycle correspondant à l’ellipsoïde d’incertitude E1 est propagée de sorte à obtenir l’erreur après l’étape 1E1 de propagation correspondant à l’ellipsoïde d’incertitude E2. Les erreurs de position sont ici supposées indépendantes afin de simplifier l’illustration, et la projection en bas à gauche de la montre l’absence de corrélation entre et (les ellipses à trait continu épais E2 et à trait continu fin E1 sont superposées en projection dans le plan ). Il ressort clairement de la que les incertitudes sur les erreurs de position après propagation sont inchangées. En revanche, une propagation entre et ou entre et est présente grâce à la corrélation entre ces grandeurs créée par la présence d’un gradient du vecteur de pesanteur. Autrement dit : De plus, dans cet exemple, est choisi égale au vecteur , c’est-à-dire selon l’axe . Les axes définis par et sont donc particuliers. En effet, la propagation sur d’une erreur de position donnée selon un azimut donné est nulle selon l’axe défini par et maximum selon l’axe défini par . La présente invention se distingue notamment de l’état de la technique par les éléments pris en compte lors de cette étape 1E1 de propagation et en particulier la prise en compte d’au moins une composante de la matrice de gradients spatiaux du vecteur de pesanteur, de préférence plusieurs composantes de la matrice de gradients spatiaux du vecteur de pesanteur, voire les neuf composantes de la matrice de gradients spatiaux du vecteur de pesanteur. Autrement dit, la matrice de propagation , et donc la matrice de transition comprend au moins un terme, de préférence une pluralité de termes dépendant chacun d’une composante de la matrice de gradients du vecteur de pesanteur. Comme illustré en introduction, l’utilisation d’une ou plusieurs composantes de la matrice de gradients du vecteur de pesanteur permet d’améliorer le modèle de propagation de l’erreur. Pour rappel, les composantes du gradient spatial du vecteur de pesanteur peuvent être représentées sous la forme d’une matrice (ici, dans le cas d’un repère dit géographique local – d’autres repères peuvent bien entendu être utilisés) : [Math. 24] Avec avec et j , est la composante du vecteur de pesanteur selon l’axe considéré, la coordonnée étant relative à la composante de la position sur l’axe sud-nord, la coordonnée étant relative à la composante de la position par rapport à l’axe est-ouest et étant relative à la composante de la position par rapport à l’axe pointant au zénith (dit repère géographique local). A partir de cette composante ou de ces composantes, la matrice de propagation (et donc de la matrice de transition ) et la matrice d’observation où est l’indice du cycle de Kalman considéré peuvent être déterminées. Plus particulièrement, au moins un composant du gradient spatial du vecteur de pesanteur est inséré dans la fonction de propagation des erreurs qui permet de calculer la matrice de propagation et, à partir de cette dernière, la matrice de transition . Autrement dit, le modèle de propagation du filtre prend en compte au moins une composante du gradient spatial du vecteur de pesanteur, de préférence issue d’un modèle enrichi. Pour rappel, la matrice de transition peut être déterminée à l’ordre 1 à l’aide de la matrice de propagation à l’aide de la relation suivante : [Math. 25] où est la matrice identité et est la matrice de propagation. Dans un mode de réalisation, la composante verticale du vecteur de pesanteur est mesurée à l’aide d’un gravimètre. Cependant, un tel moyen de mesure est encombrant et difficilement intégrable à un dispositif DI de localisation de petite taille. Aussi, dans un mode de réalisation alternatif déjà mentionné, la ou les composantes du gradient de pesanteur sont déterminées par le moyen de calcul MC à l’aide d’un modèle gravimétrique enrichi, à partir de la position estimée PE du porteur. Plus particulièrement, dans ce mode de réalisation, la détermination d’au moins une composante du gradient spatial de pesanteur comprend : la détermination, par le moyen de calcul MC, à l’aide de mesures issues du bloc de capteurs inertiels BSI et à partir d’un vecteur de pesanteur calculé à l’aide du modèle gravimétrique enrichi, de l’état de navigation du porteur incluant la position (la vitesse et l’attitude) du porteur, dite position estimée PE ; et la détermination, par le moyen de calcul MC, du gradient spatial du vecteur de pesanteur à l’aide du modèle gravimétrique enrichi et de ladite position estimée PE du porteur, par exemple à l’aide d’une fonction convertissant le modèle gravimétrique enrichi en vecteur de pesanteur et une approximation de sa matrice jacobienne, par exemple d’après la grille du vecteur de pesanteur, ou d’après les coefficients d’harmoniques sphériques. Le modèle gravimétrique enrichi peut par exemple être le modèle EGM2008 (cf. N. K. Pavlis, S. A. Holmes, S. C. Kenyon and J. K. Factor. “An Earth Gravitational Model to Degree 2160: EGM2008”. EGU General Assembly 2008, Vienna, Austria, April 13-18, 2008). De manière plus générale, comme déjà illustré auparavant, le modèle gravimétrique enrichi utilisé peut prendre différentes formes. La forme la plus simple est constituée de cartes géographiques du vecteur de pesanteur (trois composantes) et des gradients horizontaux du vecteur de pesanteur (les six composantes en et de la matrice introduite précédemment). Il n’est cependant pas nécessaire de cartographier les neuf composantes car il n’existe que cinq composantes indépendantes du gradient spatial de pesanteur. De manière alternative, le modèle de la pesanteur peut prendre une forme compressée, c’est-à-dire se présenter sous la forme d’une liste de coefficients, par exemple des coefficients d’harmoniques sphériques du potentiel de gravitation. D’autres bases telles que les harmoniques ellipsoïdales ou les ondelettes sphériques peuvent être utilisées. Pour les exploiter, il est possible de calculer les fonctions de la base à la position estimée selon les méthodes mathématiques de l’état de l’art, et de faire la combinaison linéaire adéquate pour obtenir le vecteur de pesanteur et la matrice des gradients spatiaux du vecteur de pesanteur. De nombreux modèles mondiaux du potentiel de gravitation décomposé en harmoniques sphériques existent. Par exemple le modèle EIGEN-6C4 ou le modèle XGM2019. Des explications sur la décomposition en harmoniques sphériques du potentiel de gravitation peuvent être trouvées dans le livre « Physical geodesy , second edition » – Bernhard Hofman-Wellenhof, Helmut Moritz, Springer Wien Network, 2006. Ce livre ne traite pas des gradients spatiaux du vecteur de gravitation ou de pesanteur. En revanche, cet aspect est par exemple traité dans l’article « On the computation of the gravitational potential and its first and second order derivatives » , R. Koop and D. Stelpstra – Manuscripta geodaetica – Volume 14 – pages 373-382. Dans un mode de réalisation, le dispositif DI de localisation inertielle comprend un moyen de calcul MC, un bloc de capteurs inertiels BSI et au moins un capteur externe CE, le procédé comprenant de préférence, avant l’étape de propagation, une phase d’alignement lors de laquelle le moyen de calcul MC initialise l’état de navigation et son incertitude. Dans un mode de réalisation, chaque cycle est divisé en intervalles de temps , ainsi qu’en P intervalles de temps de sorte que où et et sont des entiers positifs non nuls et où est la période du filtre. Dans ce mode de réalisation, comme illustré à la , lors de l’étape 1E1 de propagation, chaque cycle comprend, pour chaque intervalle repéré par l’indice , compris entre 1 et N, une sous-étape 1E11 de calcul de l’état de navigation (incluant la position, la vitesse et l’attitude du porteur) par propagation de l’état à partir des mesures du bloc de capteurs inertiels BSI, du modèle gravimétrique enrichi à la position de l’état , et des lois de la dynamique de Newton. De plus, lors de cette étape, une intégration par rapport au temps de la matrice de propagation est mise en œuvre, cette intégrale étant utilisée ensuite pour le calcul d’une matrice de transition élémentaire à intervalles réguliers (notés intervalle ). En effet, comme évoquée dans la suite, un calcul d’une matrice de transition élémentaire est effectué aux intervalles , avec où est la partie entière de x, initialisée à la matrice identité quand , étant un entier entre 1 et P, et complétée par l’intégration par rapport au temps de la matrice de propagation sur les intervalles de temps composant l’intervalle d’indice , où est la matrice de propagation correspondant à la fonction de propagation f linéarisée ( pouvant être simplifiée ou considérée constante sur l’intervalle ). Plus particulièrement, l’intégration de la matrice de propagation est réalisée à l’aide d’une version discrétisée de la relation [Math. 25] introduite précédemment, par exemple à l’aide d’un compteur initialisé à chaque intervalle entre 1 et N tel que soit multiple de positif ou nul. Autrement dit, pour considéré constant sur , la valeur est cumulée à la matrice identité au fur et à mesure des itérations de sorte à constituer progressivement la matrice de transition sur un intervalle . Pour cela, dans un exemple de réalisation, une lecture du vecteur gravimétrique enrichi à l‘itération est effectuée dans le modèle gravimétrique enrichi sous forme d’une matrice colonne exprimée dans un repère [q]. Puis, le vecteur gravimétrique enrichi est converti en vecteur de pesanteur enrichi dans le repère [q], le vecteur gravimétrique enrichi étant une fonctionnelle du vecteur de pesanteur enrichi. Enfin, à l’aide d’une matrice de passage , le vecteur de pesanteur enrichi est exprimé dans le repère [p]. A cette fin, dans un exemple de réalisation, le repère de mécanisation de la navigation est un repère [p] et le modèle gravimétrique est exprimé dans un repère [q]. Comme déjà mentionné, pour cette propagation la matrice de transition est obtenue d’après la matrice de propagation calculée selon l’état de l’art, et comprenant la sous matrice 3x3 de propagation qui propage le vecteur d’erreur de position exprimé dans le repère [p] sur le vecteur d’erreur de vitesse exprimé dans le repère [p], avec : où est la matrice de passage du repère [q] dans lequel est exprimé le vecteur de pesanteur au repère de mécanisation [p], où est la matrice transposée de , où est la matrice des gradients spatiaux approchée de la matrice colonne des coordonnées dans le repère [q] du vecteur de pesanteur enrichi calculé d’après le modèle gravimétrique enrichi, comprenant au moins un gradient spatial de ce vecteur de pesanteur, et où est une matrice 3x3 modélisant l’effet de la courbure de l’ellipsoïde de référence sur l’erreur de vitesse et dépendant du système de coordonnées dans lequel est exprimé le modèle gravimétrique. Le calcul de dépend de l’exploitation qui est faite du modèle gravimétrique enrichi. Si le modèle gravimétrique est exploité en repérant la position des points par leurs coordonnées cartésiennes dans le repère géocentrique [t] lié à la planète (voir figure à la fin), alors le repère [q] est égal au repère [t] et . Si le modèle gravimétrique est exploité en repérant la position des points par leurs coordonnées géographiques relativement à l’ellipsoïde de référence où L est la latitude géographique, G la longitude et l’altitude ellipsoïdale, alors le repère [q] est égal au repère géographique local [g] (Nord, Ouest, Haut), et : [Math. 27] où et sont les rayons Nord et Est de l’ellipsoïde de référence par rapport auquel sont définies les coordonnées géographiques, et où est l’altitude ellipsoïdale, avec : Où est le rayon équatorial de l’ellipsoïde de référence, et son excentricité. Dans un mode de réalisation, et peuvent être forcés à zéro dans la matrice qui est alors donné par : Dans un autre mode de réalisation, le repère géographique local peut être défini par les axes (Nord, Est, Bas) où et peuvent ou non être forcés à zéro, la nouvelle formule devenant : [Math. 30] Dans un mode de réalisation, lors de l’étape 1E1 de propagation, pour chaque intervalle d’indice multiple non nul de , à l’issue de la sous-étape 1E11, une matrice de transition étant disponible avec compris entre 1 et P, indice d’un intervalle , une sous-étape 1E12 comprend le calcul de la matrice de transition initialisée à la matrice identité au début du cycle k, obtenue progressivement à chaque nouvelle valeur de par le produit matriciel . Pour rappel, la matrice de transition pour le cycle est obtenue par le produit des matrices de transition élémentaires de tous les intervalles du cycle k. Aussi, en réalisant le produit des matrices élémentaires à chaque nouvelle valeur de , la matrice de transition associée au cycle est donnée par la valeur de obtenue au dernier intervalle du cycle considéré. Dans un mode de réalisation, lors de cette sous-étape 1E12, l’erreur de navigation et sa matrice de covariance sont également calculé à partir de l’état d’erreur de navigation de l’intervalle précédente à l’aide de la relation suivante : Où est une matrice de propagation élémentaire, étant donné par et étant donné par . Ces intervalles de temps permettent d’entretenir la corrélation qui existe entre l’erreur de position et l’erreur de vitesse pendant le cycle considéré, jusqu’au moment de l’observation effectuée lors de l’étape 1E2 de mise à jour (décrite dans la suite). L’étape 1E1 de propagation comprend également, à l’issue du dernier intervalle ( et ), une sous-étape 1E13 de calcul de la matrice d’observation à partir de l’état de navigation propagé et d’un modèle d’observation, et le calcul de propagation de l’état d’erreur et de la matrice de covariance de cet état : Dans un mode de réalisation, le dispositif DI de localisation inertielle comprend un gravimètre mesurant le module de pesanteur et (outre le calcul de la matrice de transition basée sur la matrice des gradients spatiaux approchée du vecteur de pesanteur enrichi) le calcul de la matrice d’observation en fin de l’étape 1E1 de propagation est réalisé à l’aide des gradients horizontaux de la composante verticale de pesanteur dans le repère [g] extraits du modèle gravimétrique enrichi. Dans un mode de réalisation, la matrice de transition comprend également des termes relatifs au vecteur de pesanteur (en plus du gradient), et la valeur de ce vecteur de pesanteur est déterminé à l’aide du modèle gravimétrique enrichi. Dans un mode de réalisation, afin d’améliorer la robustesse de l’estimation, un modèle d’erreur du modèle de pesanteur peut être utilisé lors de l’étape E1 de propagation. L’étape 1E2 de mise à jour Le procédé 100 selon l’invention comprend ensuite une étape 1E2 de mise à jour. L’étape 1E2 de mise à jour comprend d’abord une sous-étape 1E21 de mesure, à l’aide du capteur externe, d’une observable relative à au moins une fonction de la vitesse et/ou à la position du porteur. Il est entendu qu’une observable de la vitesse ou de la position est considérée comme observable relative à au moins une fonction de la vitesse ou de la position – de manière générale une observable est considérée comme une fonction de la vitesse ou de la position lorsque la valeur que prend cette observable est fonction de la valeur de la vitesse ou de la position. Pour cela, dans un mode de réalisation, le dispositif de localisation inertielle comprend un capteur d’altitude, c’est-à-dire un capteur permettant de déterminer la position du porteur selon l’axe (dans le repère local). Ce mode de réalisation est particulièrement avantageux dans la mesure où une centrale inertielle comprend en général un tel capteur. Dans un mode de réalisation alternatif ou complémentaire, le dispositif de localisation inertielle comprend un gravimètre, un capteur de vitesse (par exemple basé sur l’effet Doppler), une caméra, un télémètre laser, un lidar, un radar et/ou tout autre moyen de recalage. Dans un mode de réalisation, le modèle d’erreur du filtre stochastique comporte des états complétant le vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi, ces états participant à la propagation des erreurs de vitesse et étant recalés par les observations, et constituant ainsi une mesure gravimétrique intégrée à la centrale de localisation inertielle. Dans un mode de réalisation, le porteur est un bateau et le capteur externe est un capteur virtuel. En effet, dans la mesure où l’altitude est donnée par le niveau de l’eau il n’est pas nécessaire d’effectuer une mesure. En particulier, pour un porteur en mer, la mesure de l’altitude peut être considérée comme étant égale à zéro ou bien encore lue dans une carte d’altitude ellipsoïdale du géoïde. On notera que le dispositif selon l’invention peut comprendre une pluralité de capteurs externes, une observable par capteur externe étant alors mesurée lors de la sous-étape 1E21 de mesure. L’étape 1E2 de mise à jour comprend enfin une sous-étape 1E22 de mise à jour, par le moyen de calcul (MC), de l’état du filtre stochastique à partir de la mesure effectuée lors de la sous-étape 1E21 de mesure précédente de sorte à réduire au moins en partie l’état d’erreur . Cette étape est mise en œuvre de manière identique aux méthodes de l’état de la technique. Pour rappel, dans le filtre de Kalman étendu par exemple, une innovation est calculée correspondant à l’écart entre la mesure effectuée et la mesure prédite à l‘aide de la matrice d’observation et de l’état propagé. Une matrice de gains est alors calculée et l’état est corrigé par le produit du gain par l’innovation. La matrice de covariance de l’état est corrigée d’après ce gain et d’après la matrice d’observation . Il existe un gain optimal minimisant la trace de la matrice de covariance , dépendant de la matrice de covariance propagée et de la matrice d’observation . Pour chaque cycle , à la fin de l’étape de mise à jour, une correction à l’instant appelée recalage, illustrée en , détermine selon les techniques de l’état de l’art, la valeur initiale de l’état de navigation du cycle suivant d’après l’état à ce même instant combiné avec la partie de l’état d’erreur contenant les estimations d’erreurs de navigation, et détermine la valeur initiale au cycle k+1 d’états complémentaires contribuant à la qualité des estimations de l’état de navigation d’après leur valeur à la fin du cycle k combinée avec la partie de l’état d’erreur contenant les erreurs estimées sur ces états complémentaires, puis est remis à zéro, les recalages d’états et la remise à zéro d’états d’erreur étant appliqués sur la totalité ou une partie de ces états à l’instant , les composantes d’états recalés étant dites « états en boucle fermée » et les autres « états en boucle ouverte » ; Les états de navigation combinés à la partie navigation de l’état d’erreur (dont les composantes en boucle fermée sont à zéro) forment la solution de navigation à la fin du cycle k à l’instant , associée à des intervalles de confiance extraits de la sous-matrice de covariance issue de contenant la covariance des états d’erreur de navigation. L’étape 1E2 de mise à jour est illustrée à la où l’erreur après mise à jour correspond à l’ellipsoïde d’incertitude E3. Dans cet exemple, l’observation ne se fait que selon la composante en de la vitesse. Lorsque cette observation se produit, elle est prise en compte dans l’étape de mise à jour du filtre de Kalman et, si le gain optimal de Kalman est utilisé, l’ellipsoïde d’incertitude se réduit selon l’axe tout en restant inscrit dans l’ellipsoïde E2 et en étant tangent à ce dernier à l’intersection de l’ellipsoïde E2 et des points de l’espace d’état non observés de sorte à obtenir l’ellipsoïde E3. Par ailleurs, les points de l’ellipsoïde E2 appartenant au plan passant par l’origine sont statistiquement inchangés par l’observation car ils correspondent à des vecteurs d’état appartenant au plan orthogonal à l’axe d’observation. Ces points appartiennent donc également à l’ellipsoïde E3. Ces points forment une ellipse dans le plan non représentée à l’exception de la figure en haut à gauche où ils sont représentés par deux petits cercles. Le fait que l’ellipsoïde E3 soit inscrit dans l’ellipsoïde E2 se traduit par l’apparition d’une corrélation entre erreurs de position dans le plan . L’ellipse d’incertitude de position s’est donc réduite dans la direction correspondant à la direction du gradient horizontal de pesanteur. Exemple s de capteurs externes Comme déjà mentionné, le procédé selon l’invention fait appel à un capteur externe, voire plusieurs capteurs externes, pour mesurer une observable, voire plusieurs observables. En introduction de la description, les exemples ont été illustré à l’aide d’une mesure de l’altitude ou de la vitesse selon l’axe . Cependant, comme déjà évoqué, d’autre capteurs externes peuvent être utilisés dans le cadre du procédé 100 selon l’invention. Dans un mode de réalisation, le capteur externe est un télémètre laser. L’observation consiste alors en une mesure de distance du porteur par rapport à un point (connu ou inconnu). L’erreur de position produit en particulier une erreur d’accélération observable correspondant à celle du vecteur de pesanteur en projection sur l’axe de visée . Dans le procédé selon l’invention, quand la zone d’incertitude en position est comprise dans la zone de linéarité du vecteur de pesanteur, la matrice de gradients spatiaux du vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi permet de modéliser linéairement cette erreur d’accélération en fonction du vecteur d’erreur de position en trois dimensions. Cette erreur d’accélération crée une erreur de vitesse sur cet axe de projection, ainsi qu’un écart entre la distance prédite et la distance réellement mesurée. Les gains du filtre sont appliqués à cette écart (ou innovation) pour corriger le vecteur d’état (pour mémoire le vecteur d’état comporte les erreurs de position, de vitesse et d’attitudes, ainsi que le modèle d’erreurs des capteurs inertiels, et des capteurs externes) lors de la mise à jour du filtre. Ainsi, dans le procédé selon l’invention, l’erreur de position selon l’axe de visée est corrigée directement par la mesure du télémètre tandis que l’erreur sur les deux autres composantes, celles situées dans un plan orthogonal à l’axe de visée, sont corrigées d’après la matrice des gradients spatiaux du vecteur de pesanteur calculée d’après le modèle gravimétrique enrichi dans ledit plan, les gradients apportant le plus d’information étant ceux décrivant les variations de la projection du vecteur de pesanteur sur l’axe de visée en fonction des erreurs de position dans le dit plan, établissant ainsi des corrélations entre l’erreur de vitesse sur l’axe de visée et les erreurs de position dans ce plan. Lorsque l’axe de visée est vertical, l’observation est une observation d’altitude et ces gradients correspondent aux gradients horizontaux de la composante de pesanteur. C’est alors le cas évoqué en introduction. De manière plus générale, lorsque l’axe de visée est constant, comme dans le cas de l’axe vertical, deux gradients du vecteur de pesanteur peuvent suffire. En revanche, lorsque l’axe de visée est variable, il est nécessaire de disposer des neuf gradients spatiaux du vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi de façon à pouvoir observer les erreurs de position dans le plan orthogonal à l’axe de visée, quelle qu’en soit l’orientation. L’erreur de position est naturellement observable selon une combinaison linéaire de gradients spatiaux dépendant de l’axe de visée. On notera que le même principe peut être adopté avec un Lidar. De même, certaines centrales de localisation inertielle disposent d’une caméra. Le déplacement de points fixes par rapport à la Terre dans l’image donne une information sur le déplacement du porteur parallèlement au plan image. L’erreur de position produit en particulier une erreur d’accélération observable correspondant à celle du vecteur de pesanteur en projection dans le plan image . Dans le procédé selon l’invention, quand la zone d’incertitude en position est comprise dans la zone de linéarité du vecteur de pesanteur, la matrice de gradients spatiaux du vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi permet de modéliser linéairement cette accélération en fonction du vecteur d’erreur de position en trois dimensions. Cette accélération crée une erreur de vitesse en projection dans le plan image. Cela crée un écart entre la position prédite des points fixes dans l’image et celle réellement observée. Les gains du filtre sont appliqués à cet écart pour corriger le vecteur d’état lors de la mise à jour du filtre. La composante de position dans le plan image est corrigée directement par ce type d’observation. L’autre composante de position, celle située sur l’axe de visée est corrigée d’après le gradient spatial du vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi, les gradients apportant le plus d’information étant ceux décrivant les variations de la projection du vecteur de pesanteur dans le plan image en fonction de l’erreur de position sur l’axe de visée, établissant ainsi des corrélations entre les erreurs de vitesse dans le plan image et l’ erreur de position sur l’axe de visée. Mesurer des erreurs dans le plan image permet donc de corriger des erreurs sur l’axe de visée, l’information acquise dépendant de l’amplitude des gradients spatiaux. Lorsque l’axe de visée est constant, comme dans le cas de l’axe vertical, deux gradients du vecteur de pesanteur peuvent suffire. En revanche, lorsque l’axe de visée est variable, il est nécessaire de disposer des neuf gradients spatiaux. De manière alternative, lorsque l’axe de visée est asservi sur un point fixe, cette erreur de vitesse est présente dans la commande de l’axe de visée. Cette information peut être traitée comme une mesure afin de corriger l’erreur de position sur l’axe de visée au travers de la matrice des gradients spatiaux selon l’axe de visée de la projection dans le plan image du vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi. Cette méthode peut être plus précise que la première qui dépend du nombre de pixels dans l’image. Si la caméra vise des points proches, la succession des images donne également une information sur l’axe de visée sans l’aide des gradients spatiaux du vecteur de pesanteur. La corrélation entre l’erreur de vitesse et l’erreur de position donne alors une information supplémentaire et améliore encore l’estimation. Si la caméra vise des points lointains comme des étoiles (visée stellaire), très peu d’information le long de l’axe de visée est acquise par la succession des images car les étoiles sont trop lointaines, et seule la corrélation gravimétrique entre l’erreur de vitesse et l’erreur de position en permet ici l’acquisition. Bien entendu, plusieurs caméras peuvent être utilisées sur le même principe. Dans tous les cas, un gravimètre peut être ajouté aux capteurs externes afin d’améliorer la performance de la chaîne de verticale et d’estimer par exemple le biais accélérométrique de la chaîne inertielle verticale. Ce biais constitue en effet une limite de performance de l’exploitation de la corrélation d’origine gravimétrique entre l’erreur de vitesse et l’erreur de position quand seul le capteur d’altitude est disponible. En exploitant un gravimètre externe, il est possible ainsi d’abaisser la limite et améliorer la performance. Avantages associés au procédé selon l’invention Le procédé selon l’invention permet de faire bénéficier à des porteurs rapides les bénéfices multiples de la corrélation entre l’erreur de vitesse et l’erreur de position sans utiliser obligatoirement de gravimètre. Il permet non seulement de corriger les erreurs de position, mais également de corriger les erreurs de dérive gyroscopique, les erreurs de vitesse, les erreurs d’attitude, les erreurs de cap. La illustre cela en comparant les tendances de comportement de l’erreur de position et de cap d’une centrale de localisation inertielle hybridée en altitude, ce comportement dépendant aussi de la qualité du BSI, des éléments d’amortissement mécaniques et de l’environnement thermique ainsi que de la qualité du capteur d’altitude, utilisant trois types de modèle gravimétrique sur une échelle de temps de quelques heures : Le mode par défaut correspondant au modèle normal appliqué aux vecteurs de pesanteur et aux gradients spatiaux de pesanteur qui correspond à un premier procédé selon l’état de la technique (les deux graphiques les plus à gauche sur la ) ; Le mode parfois appelé « compensation » correspondant à un modèle gravimétrique enrichi pour décrire les déviations de verticale ou le vecteur de pesanteur complet, mais où les gradients spatiaux sont décrits par le modèle normal qui correspond à un deuxième procédé selon l’état de la technique (les deux graphiques au centre sur la ) ; Le mode de corrélation gravimétrique qui correspond à la mise en œuvre d’un procédé selon l’invention (les deux graphiques les plus à droite sur la ). Erreur de position horizontale d’une centrale hybridée en altitude Dans le mode par défaut, l’erreur de position horizontale comporte des oscillations de Schuler, ainsi que des erreurs liées aux erreurs gyroscopiques et à l’erreur de cap. Pour rappel, une partie des oscillations de Schuler est due au profil des déviations de verticale le long de la trajectoire, l’autre partie étant liée aux erreurs accélérométriques et gyrométriques, ainsi qu’aux erreurs d’alignement. L’observation du capteur de position ou de vitesse ou du gravimètre permet d’observer une partie de la position horizontale au travers de la corrélation entre l’erreur de vitesse et l’erreur de position. L’évolution de ladite position horizontale le long de la trajectoire contient une l’information sur la dérive gyroscopique et sur le cap et son observation partielle permet néanmoins de réduire les erreurs sur ces grandeurs dans les zones géographiques comportant du relief, ce qui peut permettre de remplacer avantageusement le signal GNSS lorsque celui-ci est indisponible. Dans le procédé selon l’invention, l’originalité consiste à réduire les erreurs sur cette grandeur d’après une mesure d’altitude (et non d’une mesure issue du GNSS) en utilisant des gradients horizontaux du vecteur de pesanteur calculé d’après un modèle gravimétrique enrichi dans la matrice de propagation . Pour être plus précis, les autres composantes du gradient participent également, dans une moindre mesure, à réduire les erreurs : les erreurs horizontales de position se propagent non seulement sur l’axe vertical, mais également sur les erreurs de vitesse horizontales au travers des autres gradients puis se propagent ensuite sur l’axe vertical. Ainsi, ils participent à l’observation des erreurs gyro et de cap. En mesurant une autre observable que l’altitude (par exemple avec un télémètre laser utilisé de façon répétée dans le temps), les autres gradients peuvent contribuer de manière plus importante à l’observation des erreurs gyroscopiques et de l’erreur de cap. Dans le mode compensé, une partie des oscillations de Schuler est atténuée car le modèle de pesanteur à haute résolution permet de réduire l’excitation de ces oscillations due au profil des déviations de verticale le long de la trajectoire. En revanche, cela ne permet pas de réduire la pente de l’erreur de position. Dans le mode de corrélation gravimétrique entre erreur de vitesse et erreur de position, l’algorithme du filtre de Kalman permet d’estimer l’erreur de position d’après la connaissance des gradients spatiaux du vecteur de pesanteur et d’après la mesure des moyens de recalage, quelle que soit la source des erreurs de position. Erreur de cap d’une centrale hybridée en altitude Dans le mode par défaut, l’erreur de cap est peu estimée. L’explication de sa faible observabilité peut être trouvée dans la littérature. Comme le montre la , Le mode de compensation n’améliore pas l’estimation du cap. En revanche, dans le mode de corrélation gravimétrique entre l’erreur de vitesse et l’erreur de position, lorsque le porteur effectue une manœuvre, la force spécifique horizontale mesurée par les accéléromètres est projetée dans le repère plateforme avec une erreur de cap produisant une erreur de vitesse horizontale et donc une erreur de position horizontale,, et comme une partie de l’erreur de position (la projection de l’erreur dans la direction du gradient spatial de pesanteur) est observée, il est possible d’améliorer l’estimation du cap en quelques heures. Pour un sous-marin ou un bateau de surface n’exploitant pas le GNSS, des mesures peuvent être effectuées régulièrement sur plusieurs jours. L’erreur de cap est l’un des contributeurs d’une oscillation de position horizontale de période 24 heures. L’observer sur cette échelle de temps longue permet de mieux la décorréler des autres erreurs. Autres erreurs Dans le mode de corrélation gravimétrique entre l’erreur de vitesse et l’erreur de position, les erreurs de dérive gyroscopique peuvent être estimées partiellement. En effet, les dérives gyroscopiques impactent fortement l’erreur de position horizontale. Leur signature dépend de la trajectoire du porteur mais comportent, au moins à l’échelle de quelques heures, une pente d’erreur de position, ainsi que des oscillations de Schuler. A plus long terme, elles produisent également une pente d’erreur de longitude et des oscillations de période de l’ordre de 24 heures. De plus, l’observation des oscillations de Schuler de la position donne également de l’observabilité aux erreurs accélérométriques. Résolution de problèmes d’observabilité Le procédé selon l’invention augmente la dimension de l’espace observable de façon plus ou moins marquée selon l’amplitude des gradients spatiaux du vecteur de pesanteur à la position estimée. Aussi, le procédé selon l’invention permet de résoudre des problèmes de fausse observabilité classiques. Dans ces problèmes où le modèle normal de pesanteur est utilisé et où l’altitude ellipsoïdale est faible, le vecteur de pesanteur est considéré de module quasi-constant et orienté selon la verticale géométrique en tout point de la Terre. Un changement d’orientation ou une translation ne change pas les projections du vecteur de pesanteur dans le repère géographique local. Dans des situations d’observation stationnaire où les axes non observables dépendent de l’état courant, un filtre de navigation de type Kalman peut être amené à créer de la fausse observabilité et devenir incohérent : ses covariances estimées diminuent de façon injustifiée et ne couvrent plus les erreurs possibles. Dans l’état de la technique, ces problèmes sont résolus dans certains cas par le filtre invariant ou le filtre OCEKF (de l’anglais Observability-Constrained Extended Kalman Filter – dans ce type de filtre, on utilise un modèle des axes non observables et l’on force le filtre à ne pas faire d’observation selon ces axes). Cependant, le filtre OCEKF possède un domaine d’utilisation plus étroit qu’un filtre de Kalman étendu. Le procédé selon l’invention permet également de résoudre certains de ces problèmes tout en conservant le filtre de Kalman qui a l’avantage de fonctionner dans un domaine d’utilisation plus large que le filtre invariant ou l’OCEKF. Un exemple d’un problème d’observabilité peut être donné dans le contexte de la localisation et cartographie simultanées (ou SLAM de l’anglais Simultaneous Localization and Mapping), lié au problème classique du suivi de points caractéristiques dans l’image afin de déterminer une estimation de la vitesse et de la position à partir de l’image en question. Des traitements appropriés permettent de sélectionner ces points et d’indiquer que leur vitesse est nulle. D’après le déplacement de ces points dans l’image, il est possible de déduire le mouvement de la caméra et donc celui du porteur. Cela peut être fait avec une centrale de localisation inertielle couplée avec la caméra (autrement dit, l’observable est mesuré à l’aide d’une caméra). Il est connu que lorsque le mouvement est faible, des problèmes de fausse observabilité peuvent se produire si le filtre stochastique est un filtre de Kalman étendu. Cela s’explique par le fait que l’observation est relative : la même scène peut être créée en tournant à la fois la caméra et les points caractéristiques sans changer les observations. Cela est vrai si le vecteur de pesanteur est orthogonal à l’ellipsoïde de référence et de module quasi-constant en tous points. Avec un procédé selon l’invention, la matrice des gradients spatiaux du vecteur de pesanteur permet de distinguer les différentes orientations possibles de la scène, et l’axe non observable ou faiblement observable avec les gradients du modèle normal peut devenir observable avec les gradients du vecteur de pesanteur calculé d’après le modèle gravimétrique enrichi tel qu’utilisé dans la présente invention. Un autre problème d’incohérence du filtre de Kalman dans l’état de la technique existe lorsque le porteur est en mouvement : le mouvement apparent des points de l’image selon l’axe de visée est nul. Dans des scénarios où la trajectoire est rectiligne uniforme, des incohérences peuvent alors être créées par le manque d’observabilité selon l’axe de visée. Une solution est apportée par le procédé selon l’invention par la modélisation des gradients spatiaux du vecteur de pesanteur d’un gravimétrique enrichi, ce dernier apportant de l’observabilité selon l’axe de visée et réduisant alors le risque d’incohérence du filtre de Kalman étendu. Dispositif de localisation inertielle Un deuxième aspect de l’invention concerne un dispositif de localisation inertielle (ou centrale de localisation inertielle) comprenant des moyens configurés pour mettre en œuvre un procédé selon l’invention. Plus particulièrement, le dispositif selon l’invention comprend un bloc de capteurs inertiels (ou BSI comprenant des accéléromètres et des gyroscopes ou des gyromètres), un moyen de calcul (par exemple une carte ASIC ou bien encore un processeur associé à une mémoire) et, de préférence, des éléments d’amortissement mécaniques sur lesquels est habituellement monté le bloc de capteurs inertiels afin de minimiser l’impact des chocs et vibrations sur la précision de localisation. Dans un mode de réalisation, le dispositif selon l’invention est configuré pour fonctionner dans ; un premier mode, dit mode d’alignement, qui correspond à l’initialisation de l’état de navigation à l’aide des mesures du bloc de capteurs inertiels et/ou d’un capteur externe ; un deuxième mode, dit mode de navigation, mode dans lequel le procédé selon l’invention est mis en œuvre par le dispositif selon l’invention. Le dispositif selon l’invention comprend également au moins un capteur externe CE, par exemple un altimètre, un télémètre laser ou bien encore une caméra. Dans un mode de réalisation, le moyen de calcul MC est associé à une mémoire, la mémoire comprenant les instructions et les données nécessaire à la mise en œuvre du procédé de réduction selon un premier aspect de l’invention. La mémoire peut notamment comprendre un ou plusieurs modèles de pesanteur et un ou plusieurs modèles de l’erreur de pesanteur. Dans un mode de réalisation, le dispositif DI de localisation inertielle selon l’invention est utilisé dans un dispositif de navigation inertielle. Procédé (100) de localisation inertielle d’un porteur à l’aide d’un filtre stochastique de type bayésien récursif mis en œuvre par un dispositif (DI) de localisation inertielle, le procédé comprenant une pluralité de cycles, chaque cycle comportant une étape (1E1) de propagation et une étape (1E2) de mise à jour, le procédé étant caractérisé en ce que, lors de l’étape (1E1) de propagation, la propagation est effectuée à l’aide d’un modèle de propagation prenant en compte au moins une composante d’une matrice des gradients spatiaux du vecteur de pesanteur fournie par un modèle gravimétrique enrichi. Procédé (100) de localisation inertielle selon la revendication précédente dans lequel le filtre bayésien récursif linéarise une loi de propagation. Procédé (100) de localisation inertielle selon l’une des revendications précédentes dans lequel le dispositif (DI) de localisation inertielle comprend un moyen de calcul (MC), un bloc de capteurs inertiels (BSI) et au moins un capteur externe (CE), et pour chaque cycle avec un entier positif non nul : lors de l’étape (1E1) de propagation, l’état de navigation en fin de cycle est calculé en fonction de l’état de navigation en début de cycle , de mesures inertielles et du vecteur de pesanteur obtenu d’après le modèle gravimétrique enrichi, l’état d’erreur est calculé en fonction de l’état d’erreur , et la matrice de covariance de l’état d’erreur est calculée en fonction de la matrice de covariance , le calcul de l’état d’erreur et de sa covariance s’effectuant à l’aide d’une matrice de transition au cycle et d’une matrice de covariance de bruit de modèle au cycle de sorte que : la matrice de transition étant déterminée à l’aide d’au moins une composante de la matrice des gradients spatiaux du vecteur de pesanteur obtenus d’après le modèle gravimétrique enrichi, une matrice d’observation étant ensuite déterminée à partir de l’état de navigation en fin de cycle et d’un modèle d’observation ; lors de l’étape (1E2) de mise à jour, l’état d’erreur après mise à jour et la matrice de covariance de l’erreur d’estimation après mise à jour sont calculés, par le moyen de calcul (MC) et à l’aide du filtre stochastique, à partir de l’état d’erreur avant mise à jour , de la matrice de covariance de l’erreur d’estimation avant mise à jour , d’une matrice de covariance du bruit de mesure , de la matrice d’observation et d’une mesure d’une observable relative à au moins une fonction de la vitesse et/ou de la position du porteur effectuée par le capteur externe au cours ou à la fin de l’étape (1E1) de propagation de sorte à réduire au moins en partie l’état d’erreur après mise à jour. Procédé (100) de localisation inertielle selon la revendication précédente dans lequel, chaque cycle est divisé en intervalles de temps , ainsi qu’en P intervalles de temps de sorte que où et et sont des entiers positifs non nuls et où est la période du filtre, l’étape (1E1) de propagation de chaque cycle comprenant, à partir d’une matrice de propagation : Pour chaque intervalle repéré par l’indice , une sous-étape (1E11) de calcul de l’état de navigation par propagation de l’état à partir des mesures du bloc de capteurs inertiels, du modèle gravimétrique enrichi à la position de l’état ; ; Pour chaque intervalle repéré par l’indice compris entre , et , étant un entier entre 1 et P désignant l’indice d’un intervalle avec où est la partie entière de x, une sous-étape de calcul d’une matrice de transition élémentaire initialisée à la matrice identité au début de l’intervalle de temps repéré par l’indice et complétée par intégration de la matrice de propagation par rapport au temps sur les intervalles de temps composant l’intervalle , la matrice obtenue par l’intégration de la matrice de propagation sur un intervalle de temps étant ajoutée à la matrice obtenue lors de l’intégration sur l’intervalle de temps précédent de sorte à constituer progressivement cette matrice de transition élémentaire , et une matrice de transition sur le cycle k initialisée à la matrice identité étant calculée progressivement à chaque nouvelle valeur de par le produit matriciel ; A l’issue du dernier intervalle , une sous-étape (1E13) de détermination de la matrice d’observation à partir de l’état de navigation propagé et d’un modèle d’observation, et de propagation de l’état d’erreur et de la matrice de covariance à l’aide de la matrice de transition associée au cycle considérée et donnée par la matrice élémentaire calculée à l’intervalle et de la matrice de bruit de modèle à l’aide des relations : Procédé selon l’une des revendications précédentes dans lequel le dispositif (DI) de localisation inertielle comprend un télémètre laser et/ou un lidar, et une mesure de la position et/ou de la vitesse du porteur est réalisée par télémétrie lors de l’étape (1E2) de mise à jour. Procédé selon l’une des revendications précédentes dans lequel le dispositif (DI) de localisation comprend une caméra et une mesure de la position et/ou de la vitesse du porteur est réalisée par la caméra lors de l’étape (1E2) de mise à jour. Procédé selon l’une des revendications précédentes dans lequel le dispositif (DI) de localisation comprend un moyen de mesure de la position verticale et une mesure de la position et/ou de la vitesse verticale du porteur est réalisée par ledit moyen de mesure lors de l’étape (1E2) de mise à jour. Dispositif (DI) de localisation inertielle comprenant un bloc de capteurs inertiels (BSI), au moins un capteur externe (CE) et des moyens configurés pour mettre en œuvre un procédé selon l’une des revendications précédentes. Programme d'ordinateur comprenant des instructions qui conduisent le dispositif (DI) selon la revendication précédente, lorsque ces instructions sont exécutées par le dispositif (DI), à mettre en œuvre le procédé (100) selon l’une des revendications 1 à 7 lorsque ledit programme est exécuté par un ordinateur. Support lisible par ordinateur, sur lequel est enregistré le programme d'ordinateur selon la revendication précédente.