Brevet FR3094082A1 — publié le 20200925

Description Titre de l'invention : Procédé d'estimation des valeurs d'harmonisation d'un magnétomètre installé dans un engin mobile, dispositif et programme d'ordinateur associés [OOO1 j La présente invention concerne le domaine de la détermination du cap magnétique d'un engin mobile à partir de mesures du champ magnétique terrestre par un magnétomètre installé dans l'engin mobile. [0002] Plus particulièrement, l'invention se rapporte à un procédé d'estimation des valeurs d'harmonisation d'un magnétomètre installé dans un engin mobile, ledit magnétomètre étant associé à son repère de référence XM, YM, ZM et l'engin mobile étant associé à son repère aéronef XA, YA, ZA, ledit procédé comprenant une étape d'acquisition d'une pluralité de valeurs de vecteurs de champ magnétique mesurées par le magnétomètre correspondant à des positionnements successifs de l'engin mobile selon des caps successifs variés. [0003] Dans l'élaboration d'un cap magnétique, la connaissance de l'orientation, dans un repère de référence, du magnétomètre tel qu'installé est nécessaire géométriquement, sinon le calcul de cap magnétique est erroné. En effet, Je calcul de cap magnétique doit être élaboré à partir de la projection dans un plan horizontal du repère géographique local des composantes de champ magnétique mesurées dans le repère magnéto-métrique. [0004] Le repère du magnétomètre (plus précisément le trièdre de mesures du capteur ma- gnétomètre) une fois installé peut ne pas être aligné avec le repère de référence, ce qui introduit des erreurs angulaires dans le calcul du cap magnétique, la dégradation s'ajoutant à la non-horizontalité de piste dans le cas de mesure au sol. [0005] Pour connaître ces défauts, dits défauts d'harmonisation, il existe une procédure de mesure lourde par visée optique sur l'engin mobile, qui impose des vérifications à plusieurs étapes lors de la fabrication des tronçons dans le cas d'un avion. La première alternative est de contraindre la précision d'installation du support du magnétomètre sur la structure de l'aéronef (contrainte lourde et pas faisable sur certaines zones). La seconde alternative est d'identifier ces harmonisations lors d'une compensation en vol complète (contrainte de procédure en vol à ce jour avec un profil prédéfini en vol). [0006] Le but de l'invention est de proposer une solution qui permette de calculer simplement, en diminuant de façon notable le temps et la complexité, les angles d'harmonisation horizontale. [0007] A cet effet, suivant un premier aspect, l'invention propose un procédé d'estimation des valeurs d'harmonisation d'un magnétomètre installé dans un engin mobile, du type 2 précité caractérisé en ce que qu'il comprend en outre une étape de détermination de l'harmonisation horizontale (hx, hy) où hx correspond à l'angle d'Euler associé à la rotation autour de XA, où hy correspond à l'angle d'Euler associé à la rotation autour de YA pour passer du repère de référence XM, YM, ZM du magnétomètre au repère avion XA, YA, ZA comprenant les étapes suivantes : [0008] - estimation d'un premier angle égal à l'angle entre XM et le grand axe d'une ellipse définie dans le plan comprenant les deux axes XM, ZM par lesdites valeurs de champ acquises ; [0009] - détermination de hy comme égale audit premier angle estimé ; [0010] - estimation d'un deuxième angle égal à l'angle entre YM et le grand axe d'une ellipse définie dans le plan comprenant les deux axes YM, ZM par lesdites valeurs de vecteurs de champ acquises ; [0011] - détermination de - hx comme égale audit deuxième angle estimé. [0012] L'invention permet ainsi de réaliser l'estimation de l'harmonisation au sol, notamment pendant la phase de calcul de la compensation des défauts magnétiques de fers durs et fers doux. [0013] Dans des modes de réalisation, le procédé d'estimation des valeurs d'harmonisation d'un magnétomètre installé dans un engin mobile suivant l'invention comporte en outre une ou plusieurs des caractéristiques suivantes : [0014] - l'estimation de l'angle parmi le premier ou le deuxième angle comprend les étapes de : [0015] - détermination des coordonnées du centre de l'ellipse par la demi-somme, sur l'axe des abscisses du plan de l'ellipse, de l'abscisse minimim et de l'abscisse maximum parmi les valeurs de champ acquises et par la demi-somme, sur l'axe des ordonnées du plan de l'ellipse, de l'ordonnée minimim et de l'ordonnée maximum parmi les valeurs de vecteurs de champ acquises ; [0016] - détermination dans le plan de l'ellipse du point, parmi les valeurs de vecteurs de champ acquises, à distance maximale du centre de l'ellipse ; [0017] - calcul de l'angle en fonction desdites coordonnées déterminées du centre de l'ellipse et des coordonnées dudit point déterminé à distance maximale du centre ; [0018] - le procédé d'estimation comprend les étapes suivantes : [0019] - calcul d'une différence entre une valeur de cap de référence et une valeur de cap résultant d'une correction, en fonction des valeurs de hx et hy déterminées, d'une mesure du champ réalisée par le magnétomètre correspondant à un positionnement de l'engin mobile selon ledit cap ; [0020] déduction de la valeur d'harmonisation verticale hz en fonction de l'égalité ci- dessous, en fonction de la différence calculée et en fonction des valeurs de hx et hy déterminées : 3 [0021] hz -7- [6W - hy. osH .strup tan(?). (-sine .siup.eoslfr + coscp. sing.))) - fax. ( - sine - tant. ces& casz,b ) "(con) . cos@ - tan(I ). (sin0 . costp . costli + ging? .sinkP)), [0022] où : [0023] I est l'inclinaison du champ magnétique terrestre audit positionnement, [0024] v, 0. p sont les angles d'attitude, respectivement de lacet, tangage et roulis de l'engin mobile audit positionnement et [0025] 511-1 représentant l'erreur de calcul de cap due aux défauts d'harmonisation ; [0026] ladite égalité donnant hz est approximée par l'une des formules suivantes : [0027] hz hz 1 - tat (1)k sine); hz 6117(roscp . rose - tan (I) (sine _rose cost + sin squIP)). [0028] dans une première phase, lesdites valeurs d'harmonisation hx, hy, hz sont dé- terminées en fonction de mesures de valeurs de vecteurs de champ magnétique mesurées par le magnétomètre pendant que l'engin mobile est au sol, lesdites étapes suivantes étant mises en oeuvre dans une phase ultérieures pour préciser les valeurs d'harmonisation, l'engin mobile étant en vol : [0029] - acquisition d'une pluralité de valeurs de vecteurs de champ magnétique mesurées par le magnétomètre correspondant à des positionnements successifs de l'engin mobile en vol selon des caps successifs variés ; [0030] - approximation de la matrice d'harmonisation telle que précisée en vol [MI,mr',,,,,,].[M 1,,'''',,1 par [1\41,,'''),] + - Id ), où Id est la matrice Identité, [Mb,'),'1 = 1 hi - hy et - hi 1 hx hy - hx 1 [0031] [0032] - lesdites valeurs d'harmonisation hxvol, hyvol, hzvol étant déterminées sur la base d'un algorithme récursif, en fonction des valeurs de vecteurs de champ magnétique acquises en vol, ledit algorithme récursif exploitant que les composantes verticales de champ magnétique obtenues après harmonisation en fonction de ladite approximation dc [Mharmosoll LI" harmmoll des valeurs vecteurs de champ magnétique sont égales à une constante ; [0033] Suivant un deuxième aspect, la présente invention propose un programme [m,,,,,' 1 hzvol - hyvol ; - hzvol 1 hxvol hyvol - hxvol 1 d'ordinateur comportant des instructions logicielles qui, lorsqu'elles sont exécutées par un ordinateur, mettent en oeuvre un procédé tel que défini ci-dessus. [0034] Suivant un troisième aspect, la présente invention propose un dispositif électronique d'estimation des valeurs d'harmonisation d'un magnétomètre installé dans un engin mobile, ledit magnétomètre étant associé à son repère de référence XM, YM, ZM et l'engin mobile étant associé à son repère avion XA, YA, ZA, ledit dispositif étant adapté pour acquérir une pluralité de valeurs de vecteurs de champ magnétique mesurées par le magnétomètre correspondant à des positionnements successifs de l'engin mobile selon des caps successifs variés ; [0035] ledit dispositif étant caractérisé en ce qu'il est adapté pour, pour déterminer l'harmonisation horizontale (hx, hy) où hx correspond à l'angle d'Euler associé à la rotation autour de XA, où hy correspond à l'angle d'Euler associé à la rotation autour de YA pour passer du repère de référence XM, YM, ZM du magnétomètre au repère avion XA, YA, ZA, estimer un premier angle égal à l'angle entre XM et le grand axe d'une ellipse définie dans le plan comprenant les deux axes XM, ZM par lesdites valeurs de champ acquises, déterminer hy comme égale audit premier angle estimé, estimer un deuxième angle égal à l'angle entre YM et le grand axe d'une ellipse définie dans le plan comprenant les deux axes YM, ZM par lesdites valeurs de vecteurs de champ acquises et déterminer - hx comme égale audit deuxième angle estimé. [0036] Ces caractéristiques et avantages de l'invention apparaîtront à la lecture de la des- cription qui va suivre, donnée uniquement à titre d'exemple, et faite en référence aux dessins annexés, sur lesquels : [0037] [fig.1] La figure 1 représente une vue d'une ensemble d'opérations mises en oeuvre dans un mode de réalisation de l'invention ; [0038] [fig.2] La figure 2 est une vue d'un aéronef dans un mode de réalisation de l'invention ; [0039] [fig.3] La figure 3 représente schématiquement un système de détermination de cap 1 dans un mode de réalisation de l'invention ; [0040] [fig.4] La figure 4 est une vue illustrant les angles d'harmonisation entre le repère du magnétomètre 2 TMAG et le repère aéronef TA ; [0041] [fig.5] La figure 5 est un organigramme d'étapes mises en oeuvre dans un mode de réalisation de l'invention ; [0042] [fig.6] La figure 6 est une vue illustrant le calcul de l'angle d'harmonisation hy dans un mode de réalisation de l'invention ; [0043] [fig.7] La figure 7 est une vue illustrant le calcul de l'angle d'harmonisation hx dans un mode de réalisation de l'invention. [0044] Considérons dans ce qui suit un aéronef 10 représenté en figure 2. Dans le mode de réalisation considéré, l'aéronef 10 comporte un système de détermination de cap 1, re- présenté en figure 3 et adapté pour déterminer un cap de l'aéronef de façon autonome par rapport à l'unité de navigation déterminant le cap primaire de l'aéronef il comprend par exemple pour cc faire un instrument secondaire d'Horizon artificiel, encore appelé « Stand-By attitude » (non représenté). [0045] Le système de détermination de cap 1 selon l'invention comprend un magnétomètre 2, une Unité de Mesure Inertielle 3, un bloc de traitement 4 et une interface homme-machine 4. [0046] Il permet de faciliter une procédure de compensation magnétique et de façon induite améliorer la qualité de la compensation, enfin pour restituer au pilote le cap autonome esti nié. [0047] Le bloc de traitement 4 est adapté pour acquérir en entrée et en temps-réel les mesures de champ magnétique délivrées par le magnétomètre 2. [0048] L'Unité de Mesures Inertielles 3, comportant de façon usuelle trois accéléromètres et trois gyromètres, détermine les angles d'attitude de l'aéronef 10 (dans d'autres modes de réalisation, les angles d'attitude pourraient être déterminés en fonction d'un autre dispositif de détermination de position et d'orientation plus précise, par exemple une unité de navigation primaire, ou un cap GPS hi-antenne, orientation fournie en temps réel en entrée du bloc de traitement 4. [0049] Différents repères sont considérés dans cc qui suit : [0050] le repère propre au magnétomètre, nommé ici TMAG, d'origine d'axes X m.Ym,Zm perpendiculaires entre eux constituant un trièdre de mesure magnéto-métrique (ou un repère de montage du magnétomètre et dans ce cas relié par une matrice d'harmonisation propre au magnétomètre calibrée en usine) ; le repère Aéronef, nommé ici TA, d'origine G point de référence de l'avion, par exemple son centre de gravité ou un point de référence mécanique de la structure, d'axes XA, YA, ZA perpendiculaires entre eux usuellement où X, est l'axe longitudinal orienté vers le nez de l'aéronef (ligne de foi), YA est parallèle aux ailes et dirigé vers l'aile droite, ZA étant lui dirigé vers le bas dans le sens croissant, comme représenté en figure 2 ; le repère Géographique Local, nommé ici TGL (encore appelé NorthEastDown) qui est lié au point courant où est l'aéronef 10 en altitude, i.e. mobile avec l'aéronef, d'origine le point de surface de la Terre à la verticale du centre de gravité G de l'aéronef, d'axes X, Y, Z perpendiculaires entre eux où le plan X, Y est un plan tangent à la surface de la Terre, X pointant vers le Nord Géographique (différent du Nord Magnétique), Z pointant vers le bas selon la verticale locale (cf figure 2) ; le repère Aéronef Horizontal (en « vol horizontal »), nommé ici THORZ, d'origine G point de référence de l'avion, par exemple son centre de gravité, 6 d'axes Xuoizz,Ynottz..ZuoRz, déduit du repère TA après correction de l'assiette de l'aéronef définie ci-après, et donc 4101 parallèle à Z du repère TGL.; [0051] Les angles d'attitude de l'aéronef sont définis par trois angles d'Euler : y (Lacet autour de ZA ou cap géographique, 0 (Tangage autour de YA) et p (Roulis, autour de XA ) comme représentés en figure 2. L'assiette de l'aéronef vs un repère horizontal THORZ est définie par les angles p et O. Ces angles d'Euler définissent la matrice de changement de repère, du repère de reference TGL vers le repère Avion à chaque instant ; les angles p et 11/ sont éléments de 1- n, + 0 est élément de ]- Tri 2, + nI 2 ]. [0052] L'attitude de l'aéronef 10 est définie par le changement de repère TGL vers le repère TA. [0053] La matrice B ((p, e, ip) des cosinus directeurs (convention aéronautique type 1-1- 3>2>1 angles de rotation Yaw IJJ> Pitch O> Roll p) est la matrice inverse suivante : [0054] .co.s(v):e:Qile cose4f.pin(e »i() sil, ,)sin(04 co:4 s neeos() Xo0) .çomy.).ço,5p)-,s4)(y.,1sinE(elsin(ç) -ços(Nip,4(e)± sin(:! s,(f) -sin(8) eo4).sin(477) .c040),cos.(4. _ . [0055] L'Unité de Mesures Inertielles ne permet pas de déterminer l'angle de Lacet y et donc complétement la matrice B, d'où le besoin d'élaborer un cap gyromagnétique autonome. Les calculs sont ainsi réalisés dans un repère horizontal THORZ à partir des angles de roulis p et de tangage 0 déterminés par L'Unité de Mesures Inertielles. [0056] Le bloc de traitement 4 est adapté pour effectuer les opérations ci-dessous décrites en référence à la figure 1, afin de déterminer la valeur du cap magnétique de l'aéronef 10 en temps réel. [0057] Par définition, le cap magnétique donne une orientation par rapport au Nord ma- gnétique terrestre. Elle est déterminée dans le plan horizontal du lieu, i.e. dans ioRy.- [0058] Le bloc de traitement 4 est ainsi adapté, dans l'exemple considéré, pour collecter les mesures de champ magnétiques effectuées sous la forme de vecteurs tridimensionnels nommés [H1J, pour ensuite effectuer une opération 200 de traitement de manière à compenser les défauts de mesures et les défauts magnétiques dans les valeurs de champ mesurées et pour délivrer le champ magnétique découlant de ces traitements sous la forme de vecteurs tridimensionnels nommés [Hcc»..w], de composantes HXcomp, HYcomp, FIZcomp. [0059] Le bloc de traitement 4 est adapté pour effectuer une opération 201 de projection du vecteur [Hcomp] ainsi délivré dans le repère horizontal THORZ correspondant à l'instant de mesure, le vecteur résultant de cette projection étant nommé [HHoRz], de composantes HXHoRz, HYHoRz, HZHORZ 7 [0060] r cos 9 sin g sin 9:1 o cos e cosesnif r H" HZ [0061] Enfin le bloc de traitement 4 est adapté pour effectuer une opération de détermination 202 du cap magnétique ipmAG selon la formule suivante : [0062] VMAG =aretan HORZ) H Y H ORZ [0063] Le cap magnétique ainsi déterminé est. ensuite utilisé par exemple comme référence typiquement d'une boucle d'hybridation gyromagnétique pour affichage ou comme référence de cap pour une fonction d'AHRS (en anglais « Attitude and Heading Reference System »), plus largement, avec différents niveaux d'hybridation dans des modes de réalisation. [0064] Les perturbateurs magnétiques en proximité du magnétomètre peuvent être modélisés comme suit : [0065] fers durs liés à la présence de matériaux ferromagnétiques aimantés ou de conducteurs électriques parcourus par des courants établis ; fers doux liés à des matériaux ferromagnétiques non aimantés qui dévient les lignes de champ conduisant à une erreur fonction de l'attitude de l'aéronef ; champs ré-induits en phase d'attitude dynamique (engendrés par courants de Foucault liés aux variations temporelles de flux magnétique dans la structure métallique de l'aéronef). 100661 Ce modèle de perturbateurs magnétique simplifié est suffisant pour des vitesses an- gulaires limitées et pour la précision de cap visé, sans tenir compte des champs ré-induits en dynamique (engendrés par courants de Foucault liés aux variations temporelles de flux magnétique dans la structure métallique de l'aéronef). 100671 Le champ mesuré perturbé dans un repère magnétomètre TMAG est modélisé comme suit, en fonction du champ magnétique théorique du lieu considéré dans le repère géographique Terrestre TGL et de l'attitude de l'aéronef 10 : 100681 fit 100691 e ' sont respectivement le champ magnétique théorique exprimé en repère TA / le champ magnétique mesuré exprimé en repère TMAG, « proche » de TA (typiquement moins de 30 mrad). 100701 ic; : paramètres de perturbation fers durs, ikl : paramètres fers doux. [0071] Un modèle général de compensation du champ magnétique mesuré dans TMAG s'en déduit : 8 [0072] = _ Itzx Àty [0074] Avec vecteur de compensation fers durs j et [0073] Ku- 2V,Ic .F Krv [0075] Kxy kxi matrice de compensation fers doux ([1C1= [kt°. [K,!'z K15 P(.1 ' Ky" JC, K 31. K, [0076] Dans le mode de réalisation particulier considéré ici, le modèle de com- pensation utilisé est volontairement réduit, défini par h, B, et [K] = K'. et utilisé pour une compensation réalisée au sol avec le champ compensé H, (proche de Ht) = K, " Hm + B,. Ceci est justifié par la non observabilité mathématique au sol et donc à la très faible qualité d'estimation de certains paramètres « verticaux » qu'il est préférable de forcer à zéro en tenant compte de l'ordre de grandeur physique des paramètres « verticaux » non diagonaux. [0077] Ainsi, en compensation au sol, une identification de 4erne ordre est appropriée avec b , by, k'=1, k'=kyx, kyy, et les approximations suivantes kn=k'=kyz=kzy=0, k'=1, bz=0. [0078] -Ken a1/4 Kcky h:RIC L 0 hmzz 100791 En approximation la matrice des fers doux est prise symétrique kxy=kyx, kyz=lczy, k'=k ' (en relation avec la déformée des lignes de champ et la physique d'un potentiel de champ dans un voisinage). 100801 Comme expliqué plus haut, concernant les défauts de mesure cette fois, il convient de déterminer, puis corriger les trois angles, dits d'harmonisation, donnant l'orientation du repère du magnétomètre TMAG (i.e. du trièdre de mesure magnétométrique du capteur), par rapport au repère Aéronef TA une fois l'équipement installé sur l'aéronef 10. 11 s'agit donc de connaître la matrice de changement de repère, matrice orthogonale dite de rotation et correspondant à la composition des trois rotations élémentaires des trois angles d'Euler en 3 dimensions. Ces trois angles, hx, hy, hz correspondant aux imprécisions dans le montage angulaire du magnétomètre sont représentés en figure 4 illustrant le passage du repère TA au repère TMAG via les rotations élémentaires 9 d'angle de valeur hz (autour axe Z), puis d'angle de valeur hy (autour axe Y), puis d'angle de valeur hx (autour axe X) avec les trièdres ainsi transformés : respectivement le trièdre [XA,YA,ZA] devient le triedre[XA',YA',ZA] après rotation de hz (autour de ZA), qui devient le trièdre [XA",YA',ZA1 après rotation autour de hy autour de YA', qui devient le trièdre [XA",YA",ZA"] après rotation de hx autour de XA" Les angles hx (angle de roulis), hy (angle de tangage) sont les angles d'harmonisation horizontale tandis que l'angle hz est l'angle d'harmonisation de lacet. Ces angles peuvent être représentés sous la forme d'une matrice dite d'harmonisation M'. [0081] Par ailleurs, étant exposé les points suivants : [0082] 1. Expression du champ théorique magnétiqueTA en repère TA à partir theo des caractéristiques du champ magnétique r dans le repère TGL (champ GL de norme HO et d'inclinaison de champ I, l'inclinaison de champ étant l'angle décrit par le champ avec sa composante horizontale, qui augmente avec la latitude sur le globe terrestre) [0083]i H.4 (hf3 = IA tho - t Te - GL "Tt' HL= 0 sin (I) [0084] 1. Projection du champ théorique magnétique en repere Aéronef Horizontal. [0085] - . TA theo HlloR B CP, , 0 )-H z...-theo = ( e W [0086] 1. Projection du champ réel magnétique en repère Horizontal tenant compte des erreurs d'harmonisation du magnétomètre monté sur Avion, définies par la matrice anti-symétrique Miu,'', : [0087] Mhart',0 -..-=1 [ hz - hy I - hz 1 hx hy - hx 1 100881 L'erreur induite par les harmonisations sur les composantes de champ s'en déduit : 100891 61-IORZ = H HORZ reel - H HORZ :herr = B ( Q, 6), (p= ()) - (41hannu - Id ) - B - H GL [0090] 1. Expression de Ferrem de cap magnétique induite par différentielle du cap : [0091] = - atan HnoRz} ) => I/ d 64) - d HHeijoRzx) - (511 oRzr [0092] Soit A 1 `1111 H 0.COS (1 ) - (SilliP - 6H HORZX CentP - ô H HORZY) [0093] 1. Soit en injectant les composantes d'erreur,_,agnoRzx et Eacilloun. fonction des harmonisations hx, hy, hz (cf point 3 ci-dessus) : HORZX dip d noRzy). 6 l H Inne7Y [0094] . cos - .0, (si osep : f ,+ SinSieP)) ± 41(0,Esee.Sire, kt hi. (-sine - tan(f); COSe. cosy [0095] Plus spécifiquement, dans le cas de l'invention, les points suivants sont démontrés par les inventeurs dans un développement d'erreur analytique et validés en simulation pour quantifier les limites d'approximations et d'harmonisation : [0096] L'harmonisation hz peut être considérée comme une erreur séparable sur le calcul de cap pour une assiette proche de l'horizontale. Cette erreur reste séparable au premier ordre après compensation au sol (le défaut hz intervient en terme anti-symétrique combiné à la matrice des fers doux et conduit après compensation au sol à des pseudos fers durs résiduels non compensés en X et Y). Cette erreur liée à hz est majeure et requiert une estimation et correction au sol. De plus pour estimer hz au sol avec une précision suffisante, hx/y doivent être assez faibles dans les limites énoncées plus loin. Les harmonisations hx et hy ont au sol un impact de second ordre après compensation au sol. Par contre leur impact devient significatif en vol en manoeuvre, d'où l'importance de corriger hx/hy pour améliorer les performances de cap magnétique en vol. On montre analytiquement dans le produit de la matrice de fers doux et de la matrice d'harmonisation, que les angles d'harmonisation peuvent être estimés au premier ordre pour le cadre énoncé ici des approximations et des limites d'harmonisation. On s'appuie sur les propriétés générales des ellipsoïdes projetées en 2D en ellipses selon deux plans : en effet le vecteur champ trace une ellipsoïde pour les points mesurés selon l'hypothèse de modèle de défauts magnétiques énoncé préalablement, le calcul de norme de champ constant localement donnant l'équation en 3 D. On énonce ici l'idée particulière que les harmonisations hx/hy contribuent en projection dans les plans respectivement IXZI/IYZI à l'angle de Faxe principal de chacune des deux ellipses et que ces harmonisations peuvent être estimées en approximation au premier ordre via l'angle de ces deux ellipses. 100971 Plus spécifiquement, dans le cas de l'invention, dans le bloc de traitement 4 est adapté pour déterminer les valeurs des angles d'harmonisation et pour corriger dans les mesures de champ magnétique mesurées par le magnétomètre, les erreurs d'harmonisation en fonction de ces valeurs. 100981 Ainsi dans un mode de réalisation, le bloc de traitement 4 comporte une mémoire 41 et un microprocesseur 42. La mémoire 41 comporte des instructions logicielles qui lorsqu'elles sont exécutées sur le microprocesseur 42, mettent en oeuvre les étapes 11 réalisées par le bloc de traitement 4 et détaillées ci-dessous en référence à la figure 5. [0099] Dans le mode de réalisation considéré, les étapes 102 et 103 ont lieu au cours de l'opération 200. [0100] Dans une phase préalable (ACQ1) 101 d'acquisition et de préparation, le champ ma- gnétique est mesuré par le magnétomètre 2 tout au long d'un tour de cap de 360' décrit par l'aéronef au sol à partir d'un cap de référence si disponible, par exemple sur une aire au sol de compensation magnétique certifiée, lieu favorable d'un point de vue magnétique. Par exemple, la rotation au sol est de 19s à 59s maximum. Ces mesures de champ magnétique (des vecteurs en trois dimensions, mesurées dans le repère du magnétomètre) sont enregistrées chacune en correspondance avec l'attitude correspondante de l'aéronef 10 telle que délivrée par l'Unité de Mesure Inertielle 3. [0101] Par exemple, un ensemble de mesures est réalisé avec acquisition à 60 Hz sur un tour au sol (plage rotation supérieure ou égale à de 360') : soit compte-tenu de rotation entre 19s et 59s, toutes les 1 minute à 5 minutes d'angle sur un tour complet. [0102] Puis un filtrage temporel des composantes de champ mesurées brutes et un même filtrage des attitudes (adapté au contexte de rotation au sol de 19s à 59s maximum) est réalisé, par exemple avec un filtre passe-bas FIR du lème ordre ou un filtre médian sur une fenêtre temporelle centrée (le cas échéant avec un filtre passe-bande en fonction des caractéristiques de bruit du magnétomètre 2). [0103] Un critère de cohérence est ensuite appliqué avec vérification d'anomalie magnétique ou de perturbations au sens large, ainsi qu'un critère de stabilité et d'horizontalité pendant la phase de mesure. (via attitude Stand-By ou via une référence externe d'attitude si disponible), pour supprimer les mesures en-onées. [0104] Le principe consiste à intégrer les variations magnétiques sous forme de critères combinés sur les composantes de champ, la norme de champ et sur l'inclinaison magnétique. L'observation bénéficie des mesures d'horizontalité et de vitesse de rotation du Stand-by attitude. Deux fenêtres temporelles sont appliquées lors de la procédure dimensionnées pour la dynamique aéronef et la classe de dérive gyrométrique. La stabilité temporelle est évaluée sur les signaux mesurés par moyenne glissante et moyenne quadratique (rms de l'anglais « root mean square ») sur ces deux horizons temporels. [0105] L'étape de filtrage et l'étape de vérification de critère sont réalisées suivant les modes de réalisation par le bloc de traitement 4 ou un autre module de traitement du système de détermination de cap 1. [0106] Les valeurs de champ magnétique (i.e. les valeurs des coordonnées des vecteurs de champ magnétique) mesurées par le magnétomètre issues de l'étape d'acquisition et préparation et les informations d'attitude associées sont utilisées par le bloc de traitement 4 pour évaluer K et i.e. le modèle de compensation, par exemple réduit 12 tel qu'évoqué ci-dessus, à l'aide des méthodes de compensation connues (par exemple : [0107] - EP0041982A1 - Crouzct - 1981 - procédé de compensation de perturbations magnétiques dans la détermination d'un cap magnétique et dispositifs pour la mise en oeuvre des procédés. IEEE-AIAA Navigation conference - Gebre Egziabher - Elkai OH - 2001 - A non-lincar two stops estimation algorithm for calibrating strapdown magnetometer. Article of Stanford University - D. Gebre Egziabhcr- OH Elkaim - ID Powels - 2002 - calibration of strapdown magnetometer in the magnetic field [0108] Ensuite, dans une étape (EST hx, hy) 102, le bloc de traitement 4 estime les valeurs d'harmonisation horizontales hx, hy : [0109] il détermine l'angle entre l'axe Y,, et le grand axe de l'ellipse E' décrite dans le plan Y,,, ZN, par les valeurs de champ magnétique mesurées par le magnétomètre issues de l'étape 101 d'acquisition et préparation, et affecte à - hx la valeur d'angle ainsi déterminée (cf. figure 7) ; il détermine l'angle entre l'axe X,, et le grand axe de l'ellipse Ex, décrite dans le plan X,,, ZN, par les valeurs de champ magnétique mesurées par le magnétomètre issues de l'étape d'acquisition et préparation, affecte à hy la valeur d'angle ainsi déterminée (cf. figure 6). [0110] Ceci dans un mode de réalisation est implémenté via un calcul algorithmique itératif de maximum/minimum de composantes en fonction de leur signe pour identifier l'angle de chacune des deux ellipses en projection dans le repère magnétométrique en 2D. Une recherche de point caractéristique de l'ellipse permet d'estimer le centre et un point de plus grande distance du centre pour calculer l'arctangente. Une vérification du second point de distance maximale est appliquée aux mesures filtrées. [0111] Par exemple dans le plan [X,,Z,,], l'angle de l'ellipse Ex, approximant hy, l'estimateur suivant est implémenté en calcul récursif (nota : l'angle de l'ellipse projetée [YZ] approxime en réalité - Hx) : [0112] . détermination des coordonnées center_elDC, center_ellZ du centre de l'ellipse par calcul des maximum et mininum des composantes de champ sur l'axe des ordonnées Z ,'', et l'axe des abscisses X,, : max(Z,,,), min(Z,), max(X,), min(X,,) : [0113] center_ellZ,(max(Z,,)+min(14)/2 center_ellX,(max(X,,)+min(Xm))/2 ; [0114] . ces coordonnées permettent d'en déduire les coordonnées du point de l'ellipse à distance maximale du centre estimé de l'ellipse : [0115] on nomme d_center= ( (((X-center_ellX).^2) + ((Z-center_ellZ).^2)) )1\0.5; [0116] on recherche le point, nommé (maxi_center) de coordonnées (Xmaxi_center, Zmaxi_center) parmi les valeurs de champ magnétique issues de l'étape 101 (ou 13 obtenu par interpolation des dites valeurs ) et qui se trouve à distance maximale du centre : [0117] maxi center=point correspondant à inax(d center); [0118] mini center=point correspondant à min(d center) ; vérification (demi grand axe et demi petit axe), i.e. déterminer l'angle de direction de l'ellipse à partir des coordonnées du point du petit axe et comparer en cohérence avec l'estimation à partir des coordonnées du point du grand axe. [0119] Calcul de l'angle direction ellipse, image du paramètre harmonisation recherché : [0120] Angle hy =atan2((Zmaxi center-center ellZ),(-Xmaxi center+center ellX)) ; [0121] Un calcul récursif similaire est par exemple en outre mené pour déterminer hx, cette fois par rapport au plan dans le plan [Y,,Z,]. [0122] Puis le bloc de traitement 4 vérifie l'hypothèse de validité aux petits angles (i.e. que hx et hy sont inférieurs à 1.5') et que la non-horizontalité de la piste est inférieure à 2'. Si non, les angles ne sont pas valides. Si oui, les valeurs de hx et hy déterminées sont considérés comme validées. [0123] Puis dans une étape (EST hz) 103, l'angle d'harmonisation hz est à présent estimé, de la manière suivante : [0124] sous-étape 103_1 : la valeur d'harmonisation IVI,'' est actualisée en fonction des valeurs de hx, hy estimées à l'étape 102 et en affectant d'abord la valeur 0 à hz : Mhar,.. = 1 [ - hz h1Z h - . -0h/1 =[ 01 01 - hy , hx v hx hy - hx 1 sous-étape 103_2 : on calcule une cucu de cap d'il) en faisant la différence entre une valeur de cap donnée Cap'f (dite de référence) correspondant à un des points de mesure de champ de l'étape 101 et la valeur de cap déterminée en fonction de la mesure du champ magnétique H = quand l'aéronef in _7-Cf pointait vers le cap Capref corrigée par la correction d'harmonisation ho- rizontale définie par hy _Ohx 1 -h: et le cas échéant la correction de [ 0 1 h compensation, i.e. on calcule dip, par la foimule : [0125] dv = Capief- Mharmo-I. K [ FÎ-- + il ] in _re f [0126] sous-étape 103_3 : on déduit la valeur d'harmonisation verticale hz de cette valeur di.p d'erreur de cap magnétique ainsi calculée, des valeurs hx, hy déterminées à l'étape 102, de l'attitude q), 19, 4) de l'aéronef correspondant à l'instant de mesure du champ magnétique ./r"-- , de l'inclinaison Tau en ref 14 point de mesure et de la formule donnée au point 5 plus haut réécrite sous la forme : [0127] hz [0128] 6W- h v.tellSe .sincp + tan (1). ( -sine .sin ço. cosep + coNcp.sinw )) - hx.(- sine-tan (1).cose.tostp )1 COSCP . CONE, - tan (I). (sine. cosq) «cosy + sin g) .sinLP)), [0129] (équation 1) [0130] La matrice d'harmonisation est actualisée à l'aide de cette valeur hz pour la poursuite des calculs des valeurs de champ compensé1'compà partir des champs mesurés H', servant ultérieurement à calculer le cap, notamment sur la base des mesures de champ magnétique H', réalisées pendant le vol de l'aéronef, à l'aide de la formule licomp= Mharmo-l- K H + III 101311 Dans un mode de réalisation, à l'étape 103 3, une des équations suivantes est utilisées à la place de l'équation I. pour déterminer hz : 101321 hz = 69) 101331 hz = 514.0 / (1 - tan (r). sine); 101341 hz = 5411 ( cose . cos° - tan (I). (sine. COSW .00Sqi Sine ssint-P) ) [0135] les deux premières équations notamment dans le cas d'angles hz limités à 1.5 deg et d'erreurs (l'horizontalité de piste inférieures à 2 deg. [0136] Dans un mode de réalisation, l'étape 103_3 est réitérée en considérant deux caps de référence Capref : cap à 00 (i.e. au Nord) et à 90°(i.e. à l'Est) (ou encore à partir de deux directions quelconques espacées de 90') ; les valeurs de hz sont calculées pour chacune de ces valeurs, puis comparées. Seulement si elles sont bien égales, la valeur correspondante de hz est retenue et la matrice d'harmonisation Mhanno est actualisée à l'aide de cette valeur pour la poursuite des calculs de champ compensé servant ultérieurement à calculer le cap, notamment pendant le vol de l'aéronef [0137] La présente invention permet donc d'obtenir des valeurs d'harmonisation de façon simple, au sol. [0138] Ceci permet ensuite dès le début du vol de l'aéronef, de calculer des valeurs de cap magnétique plus justes à partir des valeurs de champ magnétique issues des mesures réalisées par le magnétomètre embarqué, des valeurs de champ corrigées des défauts d'harmonisation. [0139] Dans un mode de réalisation, une fois les étapes 101 à 103 réalisées au sol, le bloc de traitement 4 réalise, pendant le vol de l'aéronef 10, une étape de renforcement d'évaluation de l'harmonisation. [0140] Ainsi dans une phase (ACQ2)104 d'acquisition en vol, le champ magnétique est mesuré par le magnétomètre 2. [0141] Et des valeurs plus précises d'harmonisation sont déterminées dans une étape 105 (PREC EST) comme détaillé ci-dessous. [0142] Le calcul d'harmonisation proposé ci-dessous s'appuie mathématiquement sur la sé- parabilité de la matrice d'harmonisation en somme (en tant que petits angles) d'une matrice d'harmonisation estimée au sol et d'une matrice complémentaire estimée en vol : approximation de [Mili...soll.[WIL...,..11 Par + [Miliarno,oi] -Id ), avec la matrice [MI'sol = 1 hzsol . hysol 1 , les valeurs hxsol, hysol, hzsol - hzsol 1 hxsol hysol - Insol 1 101431 étant respectivemen 1Mharmo,011= les valeurs hx, hy et hz calculées aux et 1 [ hzvol - hyvol - hzvol 1 hxvol hyvol - -hxvol I étapes 102 103. , [0144] L'estimateur est réalisé par un algorithme récursif minimisant une fonction d'erreur, par exemple un algorithme de moindre-cané récursif, à partir des mesures de champ (trois composantes Hmx/Hmy/Hmz) obtenues par le magnétomètre 2 en vol, une fois filtrées et de l'attitude de l'aéronef 10 telle que fournie par l'Unité de Mesure Inertielle 3 pour chacune de ces mesures avec application des mêmes caractéristiques de filtrage que sur les mesures magnétométriques. Selon l'état de l'art un moindre-cané consiste à résoudre résout une série d'équations linéaires sous forme matricielle ; en effet, étant posé les éléments suivants : 101451 observation construite à partir de fonctions de vecteurs de mesures sur plusieurs acquisitions (Y) = matrice d'état construite à partir de fonctions des vecteurs de mesures sur plusieurs acquisitions 1Xml * a (paramètres inconnus recherchés) ; - le meilleur estimateur â de a au sens d'erreur quadratique est : 101461 â =1Xm,* Xml -1* Xm ,* Y 101471 - A chaque itération sont calculés avec le calcul itératif d'un gain KG et de la une matrice de covariance P sans inversion de matrice ; ainsi à l'itération n 101481 ii(e) = ii:(7/- 1):-.E- :p1)-*.(17(,,,),-,v(,),*.à(!) - [0149] ?fr., = P(i-i)-- :7i):* X(iff tEP(n -1) [0150] K., (n)= p(p -1) ej17(1.,) * El ±X(n) [0151] avec fi,x(ny 4 Fol -ne goge: scalaire. , - _ [0152] La convergence de l'estimateur vers une solution â est détectée lorsque l'erreur quadratique ou la covariance sont suffisamment faibles (proche dc 16 zéro) avec possibilité d'arrêter les itérations sur de nouvelles acquisitions [0153] Dans un mode particulier dc réalisation, l'estimateur présenté ici consiste, par exemple, à identifier la composante verticale du champ magnétique constante dans le repère TGL. Cette composante est obtenue par projection des composantes de champ mesurées dans le repère TMAG, compensées du modèle dc fers ct d'harmonisation identifiés préalablement au sol, en ajoutant pour le modèle de compensation vol, trois paramètres additionnels inconnus à identifier en vol : harmonisations vol hxvol/ hyvol/hzvol (correspondant à Mharmovol). [0154] Dans ce mode dc réalisation robservabilité est favorable à l'estimation de hxvol/ hyvol, hzvol restant un défaut séparable au premier ordre comme montré dans la présente invention, aisément observable dès que l'assiette de l'aéronef varie, par exemple lors d'un premier virage en vol. [0155] [0156] où est la composante verticale du champ magnétique Hvert = constante inconnue [0157] (équation 2) [0158] HxneompensouHync.pensoyHzucompen,,,, sont les composantes x/y/z mesurées par le ma- gnétomètre, à l'instant d'échantillonnage n.Te (Te période d'échantillonnage) compensées du modèle de fers et. d'harmonisations identifiés préalablement au sol. [0159] 1-In compensai = Ksol .L + 73] et Ksol = Mhannosol K avec lit Tl et H Hyn - n tom pensai [0160] [0161] Le modèle de compensation vol ajoute les jytrametres d'harmonisation inconnus comme suit : H = Mhanno,011. H n campenvol n compensai 101621 Le Moindre-Carré est construit dans ce mode particulier en exprimant l'équation 2 sous forme matricielle avec n lignes correspondant aux n mesures échantillonnées alimentant le Moindre-Carré : [0163] Observation Y = Hve (t) live( ) [0164] avec Hve(i) - t - sinaHxn + sincpicos cp.Hyn +cosacos cp. Hzn I. compensai co ,'m il pesoi coin ,1= 1 à n, Hve(i) étant la composante verticale de champ mesurée à l'instant d'échantillonnage en volage iTe corrigée par les valeurs d'harmonisation obtenues au sol dans les étapes 102, 103. Hvert = - swe compensvo, ,÷ 'an gheos q).Hyn rompe/mot +cose.cos Hzn Com, svol Hzn - Hxn compensol" [ Hyn compensol Hzn compensol 17 [0165] - Matrice d'état [Xm] = «1(1) e2(1) e3(1) - «1(n) e2(n) e3(n) -1 101661 avec el(n)= - sine.Hyncominoi- sin Q.cose. Hut ) [0167] avec e2(n) = - stitalizit c amie ensol + cosactis (p. H feu compense!) [0168] avec c3(n) = cose.cos 9.Hyricompensa - cosasin Hzn compense?! [0169] Vecteur inconnues a = hzvol et comme indiqué plus haut , le meilleur es- hyvol hxvoi H vert timateur â de a au sens d'erreur quadratique est : [0170] â = [Xm '* Xna] '* Xm ' * Y. [0171] [0172] Ainsi : [0173] le renfort d'estimation consiste à ajouter trois inconnues complémentaires d'harmonisation, en sus des harmonisations issues de la correction sol (il ne s'agit pas d'une initialisation des paramètres d'harmonisation comme dans les méthodes existantes), la matrice d'harmonisation finale Mhanno étant la somme de Mharmosol et MHaimovol - Id). un critère d'observabilité sélectif peut être appliqué par axe selon l'attitude courante et l'écart d'attitude intégré dans le temps par exemple pour un axe vertical dans le repère TGL avec identification de la composante verticale de champ). l'estimateur présenté ici a la particularité d'utiliser les trois composantes de champ magnétique mesurées (moins sensible au bruit par rapport aux méthodes existantes utilisant la norme de champ, voire avec une seule composante comme la composante verticale); on observe ainsi sélectivement en fonction de critères d'observabilité cités plus haut la constance de chaque composante après dé-projection utilisant l'attitude courante de l'aéronef (ici, constance correspondant à différentielle temporelle nulle). [0174] Dans un autre mode de réalisation, le bloc de traitement 4 est réalisé sous forme d'un composant logique programmable, tel qu'un FPGA (de l'anglais Field Programmable Gate Array), ou encore sous forme d'un circuit intégré dédié, tel qu'un ASIC (de l'anglais Applications Specific Integrated Circuit). [0175] L'invention permet une estimation au sol des harmonisations horizontales hx, hy par analyse de la forme géométrique du champ magnétique mesuré dans les plans [Ym, et [Xm, 441, de corriger au sol en partie des défauts importants d'harmonisation (pour 18 hx, hy