Brevet FR2141984A1 — publié le 19730126

La présente invention concerne la mesure de l'intensité de la pesanteur - ou pesanteur - et plus particulièrement l'obtention de valeurs optimales de la pesanteur terrestre et des latitudes et longitudes correspondantes b partir de mesures de pesanteur, latitude et longitude dont chacune est affectée d'un bruit aléatoire. La densité des formations géologiques souterraines est une caractéristique physique qui varie beaucoup d'un type de roche à l'autre. L'intensité du champ de pesanteur engendré par une masse souterraine est, en particulier, proportionnelle à sa densité, si bien qu'une étude de la distribution des densités b partir des champs de pesanteur à la surface conduit à des informations de grande valeur en ce qui concerne la géologie souterraine, condition de procéder aux mesures physiques appropriées. Les mesures physiques comportent non seulement les mesures de pesanteur; mais aussi les déterminations des lieux géographiques où ces mesures de pesanteur sont réalisées. Par exemple, dans le cas des procédés d'exploration gravimétrique A la mer, les informations provenant d'un navire d'exploration comportent en général des indications gravimétriques, l'infornation brute provenant d'un certain nombre d'instruments de navigation et éventuellement de mesure de vitesse. Dans certains cas, l'information concernant la navigation peut etre traitée à part de l'information de pesanteur et combinée avec elle seulement lors de l'établissement de cartes. Cependant, il est nécessaire, même dans ce cas, de combiner l'information provenant de plusieurs instruments de navigation pour obtenir une bonne évaluation de la route du navire. Un gravimètre à bord d'un vaisseau qui fait route mesure une combinaison de la pesanteur et des phénomènes liés aux mouvements du navire. Les phénomènes les plus importants sont l'accélération de Coriolis et l'accélération dite centripète due aux mouvements du navire. L'amplitude de ces phénomènes est couramment dénommée "correction d'Eutvös". L'obtention de valeurs souhaitées de la pesanteur rrà l'air libre", c'est-à-dire indépendamment des mouvements du navire, et la détermination de la route du navire à partir du gravimètre et des instruments de navigation est un problème de calcul des probabilités ou, plus précisément, d'évaluation statistique. Un procédé connu d'évaluation statistique est l'application de la théorie des filtres de Weiner décrite par exemple dans l'ouvrage de Weiner, The Extrapolation, Interpolation, and Smoothing of Stationary Time Series with Engineering Applications, John Wiley & Sons, New York, 1949. Les filtres de Weiner ne sont pas utilisables dans tous les cas. Par exemple, des systèmes variant dans le temps ne peuvent etre traités convenablement par les techni ques de filtrage de Weiner actuellement connues. Par ailleurs, ces filtres de Weiner sont utilisables si l'on admet que l'ensemble est à l'état stationnaire, c'est-à-dire qu'il existe un régime permanent. Un autre procédé d'évaluation statistique, dénommé "filtrage de Kalman-Bucy" a été mis en oeuvre récemment. L'ouvrage Filtering for Stochastic Processes With Applications to Guidance, Bucy et Joseph, Tracts in Mathematics, N 23, Interscience Publishers, 1968, décrit les équations de Kalman-Bucy. On trouve aux pages 133 et suivantes de cet ouvrage l'application des filtres de Kalman-Bucy aux problèmes de navigation. On trouve à la page 134 un ensemble dans lequel un véhicule se déplace dans le champ de pesanteur terrestre. Cependant, il nty a aucune mention de détermination de valeurs de la pesanteur, de la latitude et de la longitude par les filtres de Kalman-Bucy. L'ouvrage Stochastic Processes and Filtering Theory, Andrew H Jazwinski, Academic Press, New York et Londrès, 1970, contient également une description des procédés de filtrage de Kalman-Bucy. La brochure "Kalman Filter Augmented Marine Navigation System", de Halandaris et Ozdes, de la "Satellite Positioning Corporation" décrit un système de navigation maritime. On emploie dans ce système des filtres de Ralman-Bucy pour déterminer la latitude et la longitude. Aucune des références ci-dessus ne suggère ltemFloi des filtres de Kalman-Bucy pour effectuer des déterminations optimales de pesanteur en mème temps que des évaluations optimales des latitudes et longitudes corres- pondantes. Bien que les références ci-dessus exposent en détail les principes et la théorie des filtres de Walman-Bucy, un bref résumé de cette techni que de filtrage est donnée ci-aprLs à titre d'information de base destinée à faciliter la compréhension de L'invention. Il est nécessaire, pour appliquer la théorie des filtres de Kalman-Bucy, d'établir un modèle mathématique d'un type spécial. Sous forme abstraite > il a la structure ci-après, dénommée "système dynamique excité par un bruit blanc" : Xk+1 = Ak@k + CkWk Yk = NXk + Dkvk Dans ce système d'équation Wk avec k = 1 ..., N et Vk avec k = 1 ... N sont des séries de variables aléatoires vectorielles normales indépendantes, de moyenne zéro, et des matrices covariantes d'identification n'ayant pas nécessairement la même valeur. Xk est un vecteur d'état à n x 1 dimensions à l'instant k Yk est un vecteur d'observation à m x 1 dimensions à l'instant k. Le problème de filtrage posé et résolu par Kalman et Bucy peut etre énoncé ainsi : étant donné une succession d'observations Y. avec j = 1 ... k et une distribution antérieure, de probabilité p(Xl), de l'état initial, déterminer l'évaluation optimale de l'état Xk et de sa matrice covariante terreur. Un problème plus général est le suivant : étant donné une succession d'observations Yj avec j = 1 å N et une densité de probabilités antérieure p(Xl) de la quantité X1, déterminer la valeur optimale de l'état Xk et sa matrice covariante optimale d'erreur. Dans le cas présent, N peut être inférieur, égal ou supérieur à A. Dans la terminologie courante, ces trois cas sont dénommés respectivement prévision, filtrage et régularisation. Le problème du filtrage a été récemment résolu par D.Q. Mayne et exposé dans "A Solution of the Smoothing Problem for Linear Dynamic Systems, "Automatica, Vol. 4, pp. 73-92, (1966). Selon une caractéristique importante de l'invention, des éva- luations optimales de la pesanteur terrestre ainsi que des latitudes et longitudes correspondantes sont obtenues par es mesures de pesanteur, latitude et longitude, dont chacune est affectée d'un "bruit aléatoire" de mesure. On établit selon l'invention un modèle mathématique des procédés de mesure. Ce modèle est linéarisé par rapport à une trajectoire nominale dans un espace pluridimensionnel. On établit les matrices des équations du système linéarisé et celles des équations de mesure linéarisée. Ces matrices, associées à des mesures de latitude, longitude et pesanteur sont appliquées à un appareil de traitement de l'information de Kalman-Bucy, qui fournit les évaluations optimales. Selon une autre caractéristique de l'invention, l'appareil de traitement de Kalman-Bucy -auquel on applique les matrices et les résultats de mesure- est un filtre d'égalisation. Selon une troisième caractéristique de l'invention, les matrices appliquées au filtre de Kalman-Bucy contiennent les valeurs nominales de la latitude, de la longitude et de deux composantes de vitesse. Les matrices provenant des équations de mesure contiennent des valeurs des écarts types des vitesses, des résultats de mesures de longitude et de latitude ainsi que des résultats de mesures de pesanteur. Selon une quatrième caractéristique de l'invention, les valeurs mesurées sont obtenues à l'aide d'un système de navigation à bord d'un navire. Selon une cinquième caractéristique de l'invention, les valeurs mesurées proviennent d'un ensemble aéroporté de mesure de la pesanteur, défini dans un espace à six dimensions par trois composantes de position et trois composantes de vitesse. On fait remarquer que l'utilisation nouvelle d'un système de type connu de traitement de l'information, utilisation définie par les programmes exprimant en langage symbolique les diverses opérations que comporte le procédé de l'invention, a pour effet technique de réaliser concrètement une machine nouvelle à caractère industriel certain. D'autres objets et avantages de la présente invention seront mieux compris à la lecture de la description qui va suivre d'un exemple de réali- sation et en se référant aux dessins annexés dans lesquels la figure 1 représente une opération d'exploration de la pesanteur sur un navire en mer au cours de laquelle des repères de position sont obtenus à partir d'un système de navigation maritime à satellite. Les figures 2a à 2g représentent des courbes et formes d'onde expliquant l'invention à partir d'une simulation. Selon la figure I, le bâtiment d'exploration 11 a une latitude ~ et une longitude k. La route du bâtiment 11 est décrite à l'aide de quatre quantités ou dimensions. Outre la latitude et la longitude, les deux autres quantités employées couramment sont la vitesse vers l'est VE et la vitesse vers le nord VN (comme on l'explique ci-après, l'invention est également applicable quand la trajectoire est définie par un plus grand nombre de quan tités. Par exemple, au cours d'une étude par gravimètre aéroporté, il faut six quantités ou dimensions pour décrire la trajectoire). On détermine la trajectoire du batiment par des mesures. Par exemple, comme l'indique la figure 1, on emploie un système de navigation maritime par satellites. A intervalles réguliers, on se repère par rapport aux satellites 12 et 13. Les résultats de ce repérage sont traduits par l'ordinateur de bord 14 en valeurs de latitude et de longitude. Ces valeurs sont traduites par le calcul en valeurs de vitesse vers l'est et vers le nord. Toutes ces valeurs sont enregistrées par un système d'enregistrement 15. Une autre mesure qui est effectuée par le système d'exploration est la mesure de la pesanteur "à l'air libre" par le gravimètre 16. Toutes les mesures de latitude, longitude et pesanteur sont affectées d'un bruit aléatoire. La nature statistique de ce bruit aléatoire est connue d'après la nature générale des instruments effectuant les mesures. Tous les résultats de mesure sont enregistrés à intervalles réguliers t t d'échantillonnage et l'instant d'échantillonnage est désigné par des indices. Dans le but d'exploiter ces mesures en vue de déterminations optimales, il est nécessaire de mettre en oeuvre le procède de filtrage décrit dans le tableau I ci-après. Les mesures de ~, x et g sont indiquées en 1 sur ce tableau. TABLEAU I Etablissement des matrices destinées à l'ordinateur. 1) Mesures de ~, x et g 2) Etablissement du système d'équations qui donne les résultats des mesures 3) Linéarisation de ce système par rapport à la trajectoire nominale 4) Etablissement des matrices A et C figurant dans Xk+l = AXk + CWk 5) Etablissement des matrices M et D dans l'expression Yk S MX}C + DVk 6) Application des résultats des mesures de ~, X et g à un filtre de Kalman-Bucy. La première phase de l'opération de filtrage est l'établissement du modèle mathématique de l'opération d'exécution de mesures. Le modèle mathématique du système qui réalise les mesures à bord du navire comporte un couplage entre le champ de pesanteur terrestre et le mouvement du navire. On considère tout d'abord ce mouvement du navire. Les quantités #(t) et x (t) sont la latitude et la longitude du navire à l'instant t et VE (t) et VN(t) sont respectivement les composantes N vers l'est et vers le nord de la vitesse du navire. On peut ensuite montrer que: avec Dans ces expressions a est le rayon de l'équateur terrestre et W est l'aplatissement du sphéroïde représentant approximativement la forme réelle de la terre. On admet qu'on dispose d'évaluations ou valeurs approchées des composantes devitesse désignées par VE (t) et VN (t). Ensuite, en intégrant composantes E N les équations 1 et 2 avec ces évaluations, on obtient des évaluations de position # (t) et ~ (t). La courbe (#0(t), ~ (t), VE (t), VN (t)) dans un espace à quatre dimensions est dénommée trajectoire nominale du navire. La trajectoire réelle s'écarte comme suit de la trajectoire nominale #(t) = # (t) + X1(t) (3) #(t) = ~ (t) + X2(t) VE(t) = VE (t) + X3(t) VN(t) = VN (t) + X4(t) Par exemple, si VEO et VN sont déterminées par un gyrocompas, -ou compas gyroscopique- et un loch électromagnétique, elles représentent alors le mouvement du navire par rapport à l'eau. Dans ce cas, X3 et X4 représentent les composantes vers l'est et vers le nord des courants océaniques. Etant donné que # et ~ satisfont aux équations 1 et 2 et que VE et VN sont dans les seconds membres, on en déduit que: en peut introduire alors dans le second membre de ce système les expressions # (t) + X2(t) pour Q(t), VE et X3(t) pour VE et VN et X4 pour V On reste insi avec un système de deux équations liant les quatre variables d'état X1, X2, X3 et X4. Ceséquations représentent le modèle mathématique du procédé employé pour la mesure. L'établissement de ces équations dans l'ordinateur correspond à l'opération indiquée en 2 sur le tableau I. Ces équations sont non linéaires en Xl(t), X2(t), X3(t) et X4(t). Pour pouvoir employer les techniques de filtrage de Kalman-Bucy, il est nécessaire de linéariser ces équations par rapport b une trajectoire nominale. On peut lineariser ces équations si l'on admet que les erreurs X1, X2, X3 et X4 sont petites et si l'on développe les termes non linéaires en une série de Taylor à partir de X1=X2=X3=X4=jO, c'est-à-dire par rapport à la trajectoire nominale. Ensuite, si l'on conserve seulement les termes linéaires, on obtient le système Si l'on utilise un intervalle At et remplace les dérivées par des différences finies, on obtient alors Les équations (6) représentent le système linéarisé par rapport b une trajectoire dans un espace à quatre dimensions, ce qui correspond à l'opération désignées par 3 sur le tableau I.Ces deux équations sont sous la fo ne qui convient pour l'application de la théorie des filtres de Kalman-Bucy, mais d'autres équations sont nécessaires pour les deux autres variables d'état X3 et X4 Cellesci sont en corrélation avec les courants océaniques et peuvent etre considérées comme des variables aléatoires. Le modèle le plus simple comportant suffisamment de paramètres ajustables pour être utilisable est celui d'un bruit avec corrélation exponentielle.Il est défini par les équations : Si t3 et cy sont les durées de corrélation et si a3 At et #4 #t sont les écarts types des séquences de bruit blanc W1 et W2, t3 est la durée de corrélation de la vitesse vers l'est. La durée de corrélation mesure la persistance de l'influence des vitesses antérieures sur les vitesses futures. Par exemple, si la vitesse vers l'est était de 2 noeuds il y a 2 mn à cause des courants océaniques, il est très probable que les courants océaniques ont encore une vitesse voisine de 2 noeuds dans le présent intervalle de temps, ce qui semble indiquer une durée de corrélation dépassant 2 mn. Dans la pratique, si le navire avait une vitesse de 6 noeuds, on utilisait des durées de corrélation de 10 mn. t4 est la durée de corrélation de la vitesse vers le nord. 3Lit est l'écart type de la variation aléatoire de vitesse vers l'est à l'instant tt.Ct4 tt correspond au méme écart type pour la vitesse vers le nord. Wlk est une séquence de bruit de Gauss sans corrélation dont la moyenne est zéro et la variance unité qui représente les irrégularités des courants océaniques et des vents dans la direction de l'est et W2k est la quantité correspondante pour la direction du nord. Les gravimètres de bord mesurent une combinaison de la pesanteur à l'air libre et de divers phénomènes intéressant la navigation qui seront décrits en détail brièvement. Etant donné que la variable principale est la pesanteur à l'air libre, on introduit une variable d'état pour elle : X5. On la considère par ailleurs comme affectée d'un bruit à corrélation exponentielle et définie par On peut établir maintenant les matrices nécessaires pour le filtre de Kalman-Bucy, à savoir Ak et Ck. Si l'on extrait les inscriptions des équations (6), (7) et (8), on obtient L'établissement de ces matrices dans l'ordinateur est indiqué par l'opération 4 du tableau n0 I ci-dessus. Les autres inscriptions nécessaires pour le filtre de Kalman-Bucy sont des matrices provenant d'une équation de mesure linéarisée. Si l'on considère tout d'abord le signal de sortie du gravimètre, un gravimètre idéal donnerait un signal de sortie go(t) = g(t) - E(t) (11) dans lequel g(t) représente la pesanteur et E(t) la correction d'E#tv#s. E(t) est donné avec une précision suffisante pour le travail à bord d'un navire par Ceci est expliqué plus en détail dans ouvrage de R. B. Harlan "Etstvtss Corrections for Airborne Gravemetry", Journal of Geophysical Research, Vol. 73, pp. 4675-79, (1968). Dans l'équation 12, r~ est le rayon de la terre à la latitude ~ et V~ est la vitesse tangentielle de la terre à la latitude ~.Cela donne pour l'ellipsotde de référence E = 75,03 VE cos ~ + 0,04154 (VE2 + VN2 ) (13) Si l'on remplace les variables de l'équation 13 par leurs valeurs données par l'équation 3, on obtient Si l'on forme un développement en série de Taylor on obtient E = E - 75,03 VE sin ~ X2 + (75,03 cos + 0,08308 VE ) X3 + 0,08308 VN X4 (14) + des termes d'ordre supérieur. On supprime les termes d'ordre plus élevé pour linéariser. Dans le cas présent E est la correction d'Edtvus pour la trajectoire nominale. La pesanteur g(t) est classiquement scindée en deux termes, à savoir : g(t) = gf(t) + g5(t), relation dans laquelle gs(t) est une fonction variant lentement qui tient compte des caractéristiques d'ensemble du champ de pesanteur. En général > gs est dénommé "pesanteur corrigée de la latitude".Une formule assez ancienne donnant gsest la formule internationale adoptée en 1930 gs = 978,0490 (1 + 0,0052884 sin2 ~ (15) - 0,0000059 sin 2~) cm s Avec des valeurs plus récentes des constantes indiquées dans l'ouvrage Frank D. Stacey, "Physies of the Earth", Wiley, p. 48, (1969), on obtient une meilleure formule, à savoir gs =(978,03090 + 5,18552 sin2 - 0,00570 sin2 2~) cm s-2 (16) Ces deux formules sont en fait parfaitement arbitraires.Elles sont établies à partir des premiers termes du développement du champ de pesanteur terrestre en harmoniques sphériques. Les résultats donnés par les satellites permettent de disposer de termes d'ordre plus élevé, mais on ne les emploie pas couramment. En tout cas, g, est dénommé "pesanteur corrigée de la latitude"; gf = g-gs est par ailleurs dénommé "anomalie locale de la latitude". Dans une série de Taylor, les termes du premier ordre des variables d'état sont g5 (#) = gs(# ) + gs' (# )X2 f des termes d'ordre supérieur (17) Pour linéariser les formules, les termes d'ordre le plus élevé sont supprimés. Si l'on porte 17 et 14 dans 11, on obtient alors : Ce sont les termes du premier ordre des variables d'état. On désigne par yl(t) la quantité au premier membre qui peut être déterminée à partir des signaux provenant du gravimètre et de la trajectoire nominale. L'anomalie locale de la pesanteur est désignée par X5(t), si bien que l'on a gf(t) = X5(t) (18bis) On a déjà supposé qu'il existe un bruit à corrélation exponentielle définie par l'équation 8. A partir de l'équation 18, on peut écrire dans laquelle on a inclus le terme D11 V1 pour tenir compte du "bruit" affectant les mesures. Ceci signifie que D11 est l'écart type et V1 le bruit aléatoire. L'équation 19 est ainsi sous la forme appropriée pour l'application de la théorie des filtres de Kalman-Bucy. Les deux autres équations de mesure sont basées sur le modèle le plus simple utilisable en pratique. On admet que aux "instants" k., avec j = 1... L, on dispose de "repères" désignés par #j(m) , #j(m). On admetpar ailleurs que : équation dans laquelle V2 et V3 sont des variables aléatoires de Gauss avec une moyenne nulle et une variance unité. Dans l'hypothèse où les repères sont définis avec une précision représentée par un écart type de ff mètres, on a alors relation dans laquelle a est le rayon équatorial en mètres. Ensuite, si l'on définit Y2kj et Y3kj par : on a alors Ces équations s'appliquent seulement à des instants kj pour lesquels on dispose d'un repère.Cependant, on peut les étendre à tous les instants k. On définit tout d'abord : puis on définit Y2k et Y3k arbitrairement si k # ki avec j = 1... L. Si l'on examine la solution du problème du filtrage par filtre de Kalman tel qu'il est étudié dans l'ouvrage de Mayne précité, on voit facilement que des observations "artificielles" pour k # kj ne doivent pas etre employées avec le filtre optimal. En pratique, on écrit Y2k = Y3k = 0 quand k n'est pas un instant d'observation. On arrive ensuite a la forme définitive : Y2k = M21 (k) Xlk + D22 (k) V2k (26) Y3k I M32 (k) X2k + D33 (k) V3k On peut écrire maintenant les matrices M et D sous la forme avec L'établissement des matrices 27 et 28 ci-dessus figure en 5 sur le tableau I. Avec les caractéristiques des matrices N , Ck, M et Dk figurant dans le système d'équations (9), (10),(27), (28?, les résultats des mesures peuvent etre exploités à l'aide d'un filtre de Kalman-Bucy, comme l'indiqua l'opération 6 du tableau I. Un filtre de Kalman-Bucy approprié est décrit dans l'article de Mayne cité en référence. Il est dénommé "KALQ't et il est décrit ci-après suffisamment en détail pour permettre sa mise en oeuvra. Le sigle KALQ représente aussi un programma de mise en oeuvre du filtre régulateur de Kalman pour un système constitué par un gravimètre a bord d'un navire et des repères de position approchés occasionnels. Le mécanisme effectif d'obtention des repères ne fait pas partie de ce filtre. Il est remplacé par un modèle statistique du type décrit ci-dessus. Le filtre RALQ fait intervenir un certain nombre de paramètres comme l'indiquent les équations (9), (10), (27) et (28), à savoir - - durée de corrélation de l'erreur sur la vitesse vers ltest t4 ~I durée de corrélation de l'erreur sur la vitesse vers le nord S = durée de corrélation de l'anomalie locale de la pesanteur q - écart type du bruit blanc qui détermine l'erreur sur la vitesse vers l'est - 1 écart type du bruit blanc qui détermine l'erreur de vitesse vers le nord $5 w écart type du bruit blanc qui détermine l'erreur sur l'anomalie de la pésanteur DIT = écart type de l'erreur du gravimètre cr P écart type des repères dans chaque direction. Il faut tout d'abord choisir les valeurs de ces durées de corréla tion et les écarts types. Il faut les ajuster pour "adapter" le filtre aux conditions du monde réel. Le choix des paramètres n'est pas indépendant du choix de la trajectoire nominale. Par exemple, si cette trajectoire est calculée à partir des vitesses données par le loch électromagnétique et le gyrocompas, alors C(3 et t4 représentent tout d'abord les valeurs des changements prévus des courants océaniques dans un intervalle de temps et. On utilise en général un intervalle de temps de 2 mn et, dans certaines eaux, on- peut percevoir les variations des courants pendant cette période.Par contre, si l'on considère comme nominales les vitesses calculées à partir des informations de navigation par satellite, la trajectoire nominale est beaucoup plus proche de la trajec toire vraie. Dans ce cas, on choisit des valeurs plus petites de C(3 et #4. On décrit ci-après un procédé de simulation technique pour étudier l'effet de ces paramètres sur une ambiance contrôlée. En ce qui concerne le point de départ de la simulation, on choisit arbitrairement les vitesses nominales VE (k) et VN0(k). Par exemple, on les choisit de manière à simuler une vitesse de six noeuds pour un cap de 30 La route nominale est ensuite calculée à partir des équations 1 et 2 trans formées pour des accroissements finis, en partant d'une position initiale définie par # = 0", # =45 . Les paramètres du filtre sont choisis de la manière suivante t3 =10 mn #4 =10 mn #5 =240 mn #3 = 0,25 noeud/mn (29) 4 = 0,25 noeud/mn 5 =5 milligals/mn D11 = 2 milligals = =100 mètres Des valeurs de départ ou conditions initiales, X1O, X20, X30, N40 > X50 sont nécessaires pour simuler les variables d'état. Celles-ci sont obtenues à partir de l'équation dans laquelle les X j sont choisis indépendamment de variables aléatoires de Gauss qui changent avec une valeur moyenne égale à zéro et une variance unité. Les écarts types ~j O ont les valeurs ci-après j,O # 1, 0 = 10-4 radians #2, 0 = 10-4 radians #3, 0 = 2 noeuds (29b) #4, 0 = 2 noeuds #5, 0= 100 milligaîs On emploie les équations 6, 7 et 8 pour définir le vecteur d'état pour k = 1, 2 ... 100 en utilisant des nombres de Gauss aléatoires pour les valeurs de W, lorsque les conditions initiales sont connues. Les valeurs mesurées sont ensuite données par les équations 19 et 26 en utilisant à nouveau des nombres de Gauss aléatoires pour les quantités V. On détermine ensuite la trajectoire réelle à partir des équations 3 et des signaux simulés du gra vimetre donné par l'équation 11. Ceci est représenté sur la figure 2A sur laquelle on a porté en abscisses le temps en minutes et en ordonnées la pesanteur en gals. La référence a désigne le signal de sortie régularisé du gra vimètre. La même courbe est représentée sur la figure 2B qui représente en fonction du temps les anomalies locales de la pesanteur en ordonnées. La référence b désigne la valeur estimée très proche de la valeur réelle et c la valeur nominale. Les valeurs mesurées sont ensuite introduites dans le filtre de Kalman-Bucy pour obtenir les valeurs mesurées des variables d'état. On calcule à partir de celles-ci la trajectoire estimée et la correction d'Etstvus estimée.Sur la figure 2C, on compare la correction réelle e d'Eötvos avec la correction nominale f (droite horizontale) et avec la valeur s évaluée à l'aide du filtre de Kalman-Bucy. Sur la figure 2D, on compare les valeurs réelles e et estimées X de l'anomalie locale de pesanteur. Les figures 2E, 2F, 2G représentent les résultats obtenus dans des conditions semblables pour la vitesse, le cap et la route suivie par le navire, les mêmes lettres désignant les mêmes références que ci-dessus. En ce qui concerne les formes de réalisation de l'invention, on peut avoir recours à de nombreux systèmes complexes de navigation pour obtenir les évaluations satisfaisantes de la position, de la vitesse et du cap. Des exemples de ces systèmes sont décrits dans les documents ci-après Halamandaris, H., 1969, An Integrated Positioning System, 2nd Marine Geodesy Symposium, MIS, New-Orleans, Louisiane(LUlL);Moody, A.B., 1969, High-Accuracy Offshore Navigation, 2nd Marine Geodesy Symposium, MIS, New-Orleans, Pawley, J.K., Ship Positioning With Redundant Location and Tracking Systems, Second Marine Geodesy Symposium, MIS, New-Orleans; Reinhartsen, D., 1969, The Geo Nav System, 2nd Marine Geodesy Symposium, MIS, New-Orleans. Ces ensembles comportent des systèmes électroniques de détermination de position tels que le Decca, le Lorac et le Raydist. Ces systèmes de détermination de position sont décrits dans les documents ci-après : Dean, W.N., 1966, Application of Hyperbolic Radio Systems to Marine Geodesy, lst Marine Geodesy Symposium, Columbus, Ohio (E.U.A.), pp. 115-124; Hastings, Charles E. et Comstock, Allen L., 1969, Pinpoint Positioning of Surface Vessels Beyond Line of Sight, National Marine Navigation Meeting of Institute of Navigation, San Diego, Californie (E.U.A. > . La détermination de repères de position à partir des satellites (voir figure 1) est décrite plus en détail dans les documents ci-après Kershner, R.B., (1965), Navigation by Doppler Measurement From Near-Earth Satellites, APL Technical Digest; Newton, Robert R., (1966), The Navy Navigation Satellite System, Space Research VII, North Holland Publishing Company, Amsterdam, pp. 735-763; et Stansel, T.A. (1968). The Navy Navigation Satellite System, Description and Status, Journal of Institute of Navigation, Vol. 15, NO 3. On peut obtenir des valeurs précises de la vitesse en utilisant un sonar Doppler décrit dans les documents ci-après : Farr, Harold K. et White, Donald J., 1970, A Fan Beam Doppler Sonar for Ship Navigation, prétirage, Offshore Technology Conference, Houston, Texas (E.U.A.); Marquardt Corporation, 1968, Doppler Sonar Navigation and Docking System, Publication MR 20, 440; et Turner , E.E., Thompson, B.J., et Jackson, OH., 1966, The Raytheon Acoustic Doppler Navigator, Navigation : Journal of the Institute of Navigation, Vol. 13, n0 3, pp. 210-221. On peut déterminer le cap dans de bonnes conditions à l'aide de gyrocompas. Une extension des systèmes de mesure de la pesanteur à bord d'un navire décrits précédemment est le système de mesure de la pesanteur à bord d'un aéronef. Etant donné que les aéronefs sont soumis à des accélérations verticales et horizontales plus importantes, les difficultés de détermination de la pesanteur sont plus grandes. La vitesse et le cap doivent être connus avec précision pour appliquer la correction d'Etstv#s aux résultats. Cette correction augmente avec la vitesse et dépasse 1 gal å 370 km/h. En outre, l'altitude doit être connue avec précision, étant donné que la pesanteur varie d'environ 0,3 mgal/m. Etant donné qu'il existe une intercorrélation entre la pesanteur et la navigation, le système de mesure de la pesanteur à bord d'un aéronef doit autre considéré comme une entité. On décrit ci-après un système de ce genre. MODELE D'ETAT D'UN SYSTEME DE MESURE DE LA PESANTEUR A BORD D'UN AERONEF Pour établir un modèle d'état pour un système combiné, il faut considérer un sphéroïde de référence avec un demi grand axe a et un aplatissement W. Les équations différentielles du mouvement d'un point se déplaçant à une altitude h autour de la surface dudit sphéroïde avec une vitesse V et un cap # sont relations dans lesquelles ~ et x sont respectivement les latitudes et longitude géodésiques et VE et VN sont les composantes en direction de l'est et du nord de la vitesse. VE = V sin o (32) VN = V cos e (33) Les fonctions fl et f2 ci-après contiennent la valeur de l'aplatissement sous forme de termes du premier et du second ordre 2 fl = 1 + W(2-3 sin2 #) + W (6-9 sin2 ~ + 3 sin4 ~) (34) 2 Deux équations différentielles supplémentaires mettent en corrélation l'altitude, la vitesse Vz dans le sens vertical et l'accélération verticale A z h t = Vz (V ) = A (36) zt z Le signal de sortie du gravimètre est régularisé avant enregistrement par un filtre b afin de supprimer les effets des accélérations verticales et horizontales de courte période et pour empêcher toute erreur dans le cas d'un enregistrement numérique.On doit retrancher du signal de sortie s la quantité g5 (pesanteur corrigée de la latitude), et y ajouter algébriquement la correction H d'altitude, la correction E d'E#tvtss et celle d'accélération verticale A pour obtenir l'anomalie locale de la pesanteur. z La formule internationale donnant l'intensité de la pesanteur compatible avec l'ellipsoide internationale a été adoptée en 1930. Dans le cas présent, c'est la quantité gs ou pesanteur corrigée de la latitude à savoir gs = Ge (1 + al sin2 # + a2 sin2 2~) avec : GE = 978,049 gals (37) al = 5,2884 x a2 = -5,9 x 10 6 B = latitude géodésique les quantités Ge et &alpha;1 sont déterminées par des mesures de pesanteur sur toute la surface du globe terrestre et on calcule a2 en se basant sur l'aplatissement W de l'ellipsotde international. La correction d'altitude en mgal à partir de l'elîipsotde international est donnée par l'équation 3-57 de l'ouvrage de Heiskanen, W.A. et Venig-Meinesz, F.A., 1958, The Earth and Its Gravity Field, McGraw-Hill, New York. H = (a3 + a4 sin2 #) h + a5 h2 (38) formule dans laquelle h est l'altitude en m au-dessus du niveau de la mer et &alpha;3 = 3,0877 x 10-1 &alpha;4 = -4,4 x 10-4 (38a) a5 = -7,3 x 10-8 La correction d'Eötvus en mgals est donnée par formule dans laquelle h est l'altitude en m, VE, VN sont les composantes est, ouest et sud-nord de la vitesse en m/s &alpha;6 = 7,2921.10-5 s-1, c'est-à-dire la vitesse de rotation de la terre basée sur le jour sidéral de 86 164 s. sont les vitesses au sol sont les vitesses à l'altitude h sont les vitesses au sol sont les vitesses à l'altitude h La valeur vraie de l'anomalie locale de la pesanteur gf est donnée par : gf = g - (g -E-H+Az) (40 > En ce qui concerne le modèle linéarisé, les équations du système sont non linéaires et doivent être linéarisées pour être adaptées aux formules de Kalman. Le modèle d'état est basé sur un vecteur X d'état à huit dimensions défini par les relations 41 ci-après. X1 .= - # , erreur de latitude X2 = # - # , erreur de longitude (41). X3 = V - VN , erreur sur la composante sud-nord de la vitesse o X4 = VE - VE , erreur sur la composante ouest-est de la vitesse X5 = g, anomalie locale de la pesanteur (41 > X6 = h-h , erreur sur l'altitude X7 s Vz-Vz , erreur sur la composante verticale de la vitesse X8 IA - A , erreur d'accélération verticale. z z Chaque composante, sauf X5 est la différence entre la valeur vraie et la valeur nominale (avec l'exposant zéro) du paramètre correspondant. Le vecteur Y exprimant les résultats de mesure est le vecteur d'erreur à quatre dimensions ci-après Y Y1 = g + #1 (42) Y2 = # - # + #2 Y3 = # - # + #3 Y4 = h - h + #4 formules dans lesquelles # i avec i de 1 à 4, représentent des bruits blancs gaussiens a moyenne zéro et sans corrélation mutuelle. Pour linéariser les équations 30 et 31, on les développe en séries de Taylor par rapport à la trajectoire nominale et on supprime les termes apres les dérivées partielles d'ordre 1 et on obtient équations dans lesquelles sont les vitesses au sol sont les vitesses à l'altitude h Les équations différentielles d'erreur correspondant à 36 sont les suivantes x6 = X7 , X7' = X8 (45) Les erreurs sur les deux composantes horizontales de la vitesse, l'anomalie locale de la pesanteur et l'accélération verticale sont considérées comme étant des bruits avec corrélation exponentielle et données par les équations différentielles ci-après Les termes correctifs de l'équation 40 sont linéarisés par développement en série de Taylor et limités aux termes du premier ordre, Correction de latitude s g50 + X1 gs (50) gs' = 2Ge (&alpha;1 sin # cos ~ - 2&alpha;;2 sin 2# cos 2# (51) Correction d'altitude Correction d'Ebtvus L'équation (40) devient alors Si l'on remplace gf par X5 dans l'équation 60, on obtient L'équation donnant la valeur de la pesanteur est Y1 =g - (gs - E - H + Az ) + #1 (62) dans laquelle #1 est l'erreur de mesure surla pesanteur. En portant 61 dans 62 on obtient Medèles pour accroissement finis Les équations différentielles linéaires 43 à 49 sont remplacées pour les accroissement finis par les approximations habituelles Les équations aux accroissements finis d'état et de mesure pour le modèle mathématique du procédé de génération de l'information sont X(k+1) = A(k) X(k) + C(k) # (k) (65) Y(k) = M(k) X(kJ + D(k) # #(k) (66) Les vecteurs ss (k) et a (k) sont des séquences de bruit blanc avec une moyenne nulle et des variances unité, correspondant aux vecteurs 8 et représentant des phénomènes aléatoires dans les équations 46 - 9 et 42. Les éléments différents de zéro des matrices A, C, M et D variant dans le temps sont énumères ci après, les indices k étant supprimés sur tous les termes variant dans le temps. Matrice de transition A Matrice de bruit à l'entrée C Les variances des séquences de bruit blanc à l'entrée avec i = 3, 4, 5 et 8 sont égales à (#i# #)2, on a par conséquent Matrice de sortie M C15 = C18 1 C21 = 2 C32 i3 C46 = avec si la composante de rang i de Y est mesurée à l'instant k dans le cas contraire Matrice de mesure D Cette matrice est une matrice diagonale dont les éléments Dii sont les écarts types des séquences de bruit de mesure \ , avec i = 1, 2, 3, 4. Evaluation du vecteur d'état Le modèle linéaire discontinu qui a été établi est a la base de la réalisation d'un appareil Kalman de traitement de l'information. On fait, cela étant, passer les informations discontinues å travers un filtre de Kalman de traitement en temps réel ou un dispositif de lissage -ou de régularisation- de Kalman au cours d'un traitement postérieur au vol et on obtient ainsi une évaluation optimale du vecteur d'état. Les évaluations optimales de l'anomalie locale de pesanteur et de la trajectoire de navigation (position, vitesse et altitude) sont obtenues par addition algébrique des évaluations du vecteur d'état à la trajectoire nominale, comme on l'indique dans les équations ci-après, où les évaluations sont caractérisées par un (^). # = ~ + X1 (70) X = # + X2 VN = VN + N3 V E = V0 + N4 E g = X5 h = h + X6 Vz = vu + X7 A Z A X8 Les corrections estimées de pesanteur peuvent être déterminées en. portant les évaluations de latitude, de vitesse et d'altitude, dans les équations non linéaires 37 à 3. Le filtre (ou le dispositif de lissage) de Kalman est optimal en ce sens que si l'on considère l'ensemble de la catégorie des filtres linéaires (dispositif de lissage) auquel il appartient, il réduit au minimum la trace de la matrice covariante d'erreur. (La matrice covariante d'erreur est définie T par E(x-x).(x-x) , oh E représente l'espérance mathématique et T la matrice transposée. La trace est la somme des éléments diagonaux). Si le système de génération de l'information est linéaire et si les bruits aléatoires associés sont gaussiens, alors le filtre (dispositif de lissage) de Kalman est optimal pour toutes les catégories de filtres ou de dispositifs de lissage.Etant donné qu'on applique au dispositif de lissage des données futures disponibles, alors que le filtre ne doit fonctionner qu'à partir des informations passées et présentes, c'est le dispositif de lissage qui a le meilleur rendement, autrement dit la trace de sa matrice covariante d'erreur est inférieure à la trace de la matrice covariante d'erreur du filtre, quel que soit l'instant d'échantillonnage. La trace de la matrice covariante d'erreur est une moyenne d'ensemble théorique qui est calculée. et mise en oeuvre à l'aide d'un appareil de traitement de Kalman. Etant donné qu'elle représente une moyenne d'ensemble, elle dépend de la structure du modèle mathématique du procédé de génération de l'information et non de l'information obtenue par ce procédé. Bien entendu, diverses modifications peuvent être apportées par l'homme de l'art aux dispositifs ou procédés qui viennent d'être décrits uniquement à titre d'exemples non limitatifs, sans sortir du cadre de l'invention. REVENDICATIONS 1. Système de traitement de l'information utilisant une calcula trice automatique caractérisé en ce qu'il est programmé selon un programme particulier permettant,par exemple, d'obtenir des déterminations optimales de l'intensité de la pesanteur terrestre ainsi que de la latitude et de la longitude correspondantes, à partir de mesures entachées de bruit qui comportent une mesure de pesanteur, ce programme particulier étant tel que défini ci dessous 1. Etablissement d'un modèle mathématique dudit système et du procédé associé qui permet d'effectuer lesdites mesures; 2. Linéarisation-dudit système par rapport à une trajectoire nominale dans un espace pluridimensionnel; 3.Etablissement des matrices A et C faisant partie de l'équation du système linéarisé X = AN + C/u dans laquelle X est le vecteur d'état, lu est le vecteur bruit aléatoire d'entrée et A et C sont fonction de ladite trajectoire nominale et des opéra tions aléatoires qui décrivent ledit système; 4. Etablissement des matrices M et D figurant dans l'équation de mesure I inéarisée: Y = MX + Dv dans laquelle Y est le vecteur résultat de mesures, v est le vecteur bruit aléatoire de la mesure et M et D sont liés à la trajectoire nominale, aux caractéristiques des instruments de mesure et aux procédés aléatoires qui décrivent lesdites mesures; et 5. Application desdits résultats de mesure à un appareil de traitement de Kalman-Bucy défini par les matrices A, C, M et D. 2. Système selon la revendication 1, caractérisé en ce que ledit appareil de traitement de Kalman-Bucy est un dispositif de lissage. 3. Système selon la revendication 1, caractérisé en ce que ltopéra- tion de linéarisation comprend :le développement dudit système d'équations en série' due Taylor et la suppression des termes d'ordre élevé. 4. Système selon la revendication 1, caractérisé en ce que lesdites mesures comprennent des mesures de position qui sont effectuées à l'aide d'un ensemble d'exploration à bord d'un navire et en ce que l'opération de linéari sation comprend la linéarisation dudit modèle par rapport à une trajectoire nominale dans un espace quadridimensionnel qui comprend deux composantes ou coordonnées de position et deux composantes de vitesse. 5. Système selon la revendication 4, caractérise en ce que la matrice A a la forme ci-après dans laquelle tt est la durée de l'intervalle de temps d'échantillonnage; l'indice k désigne un intervalle de temps d'échantillonnage et l'indice k41 désigne l'intervalle Lt de temps d'échantillonnage qui lui succède; de la vitesse sudwnord et 575 est la durée de corrélation pour le champ de pesanteur terrestre. 6. Système selon la revendication 5, caractérisé en ce que la matrice C est représentée par matrice dans laquelle d, #4 et Cr sont, respectivement, les variations 3 4 prévues, par intervalle d'échantillonnage, des composantes ouest-est, sud-notd de la vitesse et de la pesanteur. 7. Système selon la revendication 6, caractérisé en ce que la matrice M est représentée par : matrice dans laquelle t 2 est égal à (gs' + 75 V sin ~t); gs' est la dérivée par rapport a la latitude de la pesanteur corrigée de la latitude, 3 est égal b -(75 cos #K + 0,08308 VE ) et E 4 à-0,08308 VN . V 8. Système selon la revendication 7, caractérisé en ce que la matrice D est représentée par matrice dans laquelle D11, D22 et D33 sont les écarts types des instruments utilisés pour les mesures de pesanteur, de longitude et de latitude. 9. Système selon la revendication I, caractérisé en ce que lesdits résultats de mesure sont des résultats de mesure de position qui sont obtenus à l'aide d'un système d'exploration aéroportée et en ce que l'opération de linéarisation réalisée par ce système comporte la linéarisation dudit modèle par rapport à une trajectoire nominale dans un espace à six dimensions qui comprend trois composantes ou coordonnées de position et trois composantes de vitesse. 10. Système selon la revendication 9, caractérisé en ce que les éléments différents de zéro de la matrice A sont et chacune des quantités figurant ci-dessus a été définie antérieurement dans le texte. 11. Système selon la revendication 10, caractérisé en ce que les éléments différents de zéro de la matrice B sont les suivants relations dans lesquelles #3, #4, #5 et #6 sont les variations prévues, par intervalle d'échantillonnage, de la composante ouest-est de la vitesse, de la composante sud-nord de la vitesse, de la pesanteur et de l'accélération verticale. 12. Système selon la revendication 11, caractérisé en ce que les éléments différents de zéro de la matrice C sont les suivants avec si la composante de rang i de Y est mesurée à l'instant k dans le cas contraire et les termes ci-dessus ont été définis antérieurement, dans le mémoire descriptif. 13. Système selon la revendication 12, caractérisé en ce que la matrice D est une matrice diagonale dont les éléments Dii sont les écarts types t i des séquences de bruit des résultats de mesure avec i = 1, 2, 3 et 4.