Une fois toutes ces opérations faites, on peut traiter les pseudo-distances. Les parties en *italiques grasses* sont des sous-parties du programme que nous n'étudierons pas en détail.

## LE PROGRAMME PSEUDO

<figure><img src="image_52.png" /></figure>

Le programme PSEUDO fait appel à un certain nombre de procédures que nous allons étudier maintenant

**CARGEO:** est la procédure qui permet de passer des coordonnées cartésiennes aux coordonnées géographiques.

**GEOCAR:** est la procédure qui permet de passer des coordonnées géographiques aux coordonnées cartésiennes.

**TRAJOND:** est la procédure qui permet de faire la correction du temps de parcours de l'onde due au mouvement de rotation de la Terre durant ce trajet.
**ANGLES:** est la procédure qui calcule l'élévation et l'azimut d'un satellite vu du point stationné.

## PROCEDURES DE CORRECTIONS IONOSPHERIQUE ET TROPOSHERIQUE

**CORIONO:** est la procédure qui calcule la correction ionosphérique à apporter aux mesures brutes. Cette correction est basée sur un modèle standard d'évolution de l'ionosphère à l'échelle de l'année.

**CORTROPO:** est la procédure qui calcule la correction troposphérique à apporter aux mesures brutes. Il s'agit d'une correction standard faisant intervenir uniquement l'altitude du point.

## PROCEDURES DE SAUVEGARDES

**SAVEHORL:** est la procédure de sauvegarde des coefficients d'horloge.

**SAVECOOR:** est la procédure de sauvegarde des coordonnées.

La procédure **DIRECT** calcule la distance entre le point stationné et le satellite.
La procédure CHOL est la procédure d'inversion d'une matrice par la méthode
de CHOLEWSKI

CHOL

<figure><img src="image_53.png" /></figure>
Nous allons maintenant voir la dernière partie du programme: le traitement
de la phase.

Nous avons en notre possession un fichier de données d'observations dans le
directory "OBSERV" ainsi qu'un fichier d'éphémérides.

Dans le directory courant nous avons un fichier "trait.dat", un fichier
"dialogue.dat" et un fichier de coordonnées issues du calcul de pseudo-distance.
<figure><img src="image_54.png" /></figure>
IMPLANTATION D' UN RESEAU GEODESIQUE
AUTOUR D'UN SYSTEME DE RIFTS ACTIFS

ANNEXE 7

Le logiciel de compensation TRIANGU
## INTRODUCTION:

La Géodésie fait appel à des mesures d'éléments géométriques (distances, angles,...), réalisées avec une très grande précision ($10^{-6}$) sur une Terre physique, le Géoïde, dont la forme est proche d'un ellipsoïde de révolution. En aucun cas, pour des réseaux de dimension supérieure au km, on ne peut négliger la courbure de la terre et assimiler le Géoïde à une surface plane. On aura donc recours à un modèle ellipsoïdique de Terre sur lequel seront définis, un système de coordonnées (coordonnées géographiques: *Longitude, Latitude*) et sur lequel il faudra réduire les mesures. Cependant pour éviter des calculs trop complexes et se placer dans un système de représentation plus familier et plan, la plupart des calculs, notamment ceux d'ajustement de coordonnées seront effectués après transformation dans un système de coordonnées rectangulaires planes. On parlera alors de réduction à une *projection conforme* (qui conserve les angles). Cette transformation se fera au détriment des données mesurées, au niveau des valeurs projetées sur l'ellipsoïde, puis à la projection par une adaptation rigoureuse d'échelle et de corrections sur les angles.

Nous utilisons ici la projection UTM (Universal Transverse Mercator), utilisable dans le monde entier, mais qui a cependant l'inconvénient de conduire à de fortes réduction d'échelle.

Les altitudes des différents points du réseau ont été ajustées séparément. Le calcul des coordonnées se fera uniquement sur les coordonnées horizontales.

## METHODE D'AJUSTEMENT

On cherche à obtenir la précision maximum sur les mesures et sur les coordonnées des points qui en résultent. Pour obtenir ce résultat et estimer correctement le degré de précision obtenu on sera donc amené à effectuer des mesures surabondantes:

- pour chaque élément mesuré sur le terrain on répétera les observations afin d'estimer l'écart type sur la mesure, en estimant que les erreurs accidentelles suivent une loi gaussienne.

- on estimera du mieux possible les erreurs systématiques des instruments ou des méthodes de mesure (après correction d'éléments perturbateurs, par exemple: l'indice de l'air pour les mesures de distances).

- les mesures sur le réseau auront une configuration géométrique telle que les coordonnées d'un point seront obtenues à partir de plusieurs autres points voisins de manière surabondante.
Cette surabondance de données mesurées par rapport au nombre de coordonnées de points à estimer conduira à un ajustement de ces coordonnées qui minimise les erreurs d'observation en utilisant une méthode de moindres carrés.

Pour effectuer cet ajustement on a choisi d'utiliser la **méthode d'inversion** développée par Tarantola et Valette (1982a,b) qui est une généralisation des méthodes de moindres carrés de résolution des problèmes inverses.

**la méthode d'inversion de Tarantola-Valette :**

Tout système physique est déterminé par un certain nombre d'éléments se subdivisant en **paramètres** représentatifs du problème et en **données** mesurées. L'objectif du problème inverse est de déterminer les valeurs d'un nombre fini de paramètres décrivant le système physique à travers un certain nombre de fonctions (le modèle physique) liant ces paramètres aux données. Tarantola et Valette (1982b) considèrent le problème comme une combinaison d'informations où ils distinguent : l'information expérimentale sur les données, une information à priori sur les paramètres, et une information théorique sur la réalité physique, c'est-à-dire le modèle choisi reliant entre elles toutes les variables du problème. Ils font l'hypothèse que les erreurs d'observations, de modélisation et sur les informations à priori sur les paramètres suivent une loi gaussienne. Ils recherchent alors une solution au sens des moindres carrés qui aboutisse à une solution stable et montrent que le problème peut être résolu par un algorithme itératif permettant de trouver les valeurs des paramètres qui maximisent la densité de probabilité.

Dans le problème de géodésie traité, les **paramètres** sont les coordonnées X et Y des points du réseau et les **données** sont les distances et les angles mesurés.

Toute l'information à priori sur les paramètres est exprimée sous la forme d'un vecteur $\mathbf{p}_j$ et d'une matrice de covariance $C_p$ calculée à partir des écarts à priori $Ep(j)$, valeurs maximales des variations tolérées des paramètres. Ici, les valeurs a priori des paramètres sont des coordonnées approchées.

Les données sont les éléments mesurés, et sont exprimées de la même façon sous la forme d'un vecteur $\mathbf{d}_o$ ; chaque donnée est caractérisée par un intervalle de confiance $E_{D(i)}$ permettant de définir une matrice de covariance $C_d$.
Le modèle théorique ou **Problème direct** est l'ensemble des algorithmes qui permet, connaissant des valeurs initiales des paramètres, de déterminer des **Valeurs Calculées** des variables mesurées.

$$d_{cal} = g(p)$$

Dans ces conditions le problème inverse peut être résolu par un algorithme itératif permettant de trouver le point [P] correspondant à la densité de probabilité maximum (Tarantola et Valette, 1982b).

$$p_{k+1} = p_k + (G_k^T C_d^{-1} G_k + C_p^{-1})^{-1} \cdot (G_k^T C_d^{-1} \cdot (d_0 - g(p_k)) - C_p^{-1} \cdot (p_k - p_0))$$

où

$G_k$ est la matrice des dérivées partielles des données calculées par rapport aux paramètres.

$$G_k^{j\alpha} = \left( \frac{\partial g_j}{\partial p^\alpha} \right) p = p_k$$

L'estimation des erreurs sur les coordonnées finales s'obtient par la matrice de covariance à postériori des paramètres:

$$C_k' = (G_k^T C_d^{-1} G_k + C_p^{-1})^{-1}$$

Le calcul des valeurs propres de $C_k'$ permet alors d'obtenir les valeurs principales et les orientations des axes des ellipsoïdes d'erreur sur chaque point.

(voir le chapitre :Mise en oeuvre de la procédure d'inversion et l'organigramme correspondant.)

## COORDONNEES GEOGRAPHIQUES SUR L'ELLIPSOIDE

Un point M à la surface de la Terre peut être repéré :

- soit par ses coordonnées cartésiennes géocentriques (X,Y,Z) dans un trièdre de référence dont le centre est voisin du centre de masse O de la Terre
- soit par ses coordonnées ellipsoïdiques, définies au point m, projection du point M selon la normale Mm à un ellipsoïde de révolution centré en O :

$λ$ = longitude du point m

$Φ$ = latitude du point m

$H$ = mM = altitude ellipsoïdique du point M

le couple ($λ,Φ$) est appelé aussi coordonnées géographiques sur l'ellipsoïde considéré.

On donne ci-dessous pour mémoire les éléments de passage entre les deux systèmes :

Le plan XOZ est le plan méridien origine

$$X = N \cdot \cos\Phi \cdot \cos\lambda$$

$$Y = N \cdot \cos\Phi \cdot \sin\lambda$$

$$Z = N \cdot (1 - e^2) \cdot \sin\Phi$$

où N est la grande normale

$$N = \frac{a}{\sqrt{1 - e^2 \sin^2 \Phi}}$$

<figure><img src="image_55.png" /><figcaption>Figure. 1</figcaption></figure>

L'ellipsoïde de référence est laissé au choix de l'utilisateur qui rentre ses paramètres dans le fichier de projection. Le plus souvent on utilise l'Ellipsoïde International.

**Remarque** : Pour les réseaux d'étude des déformations de l'Ecorce Terrestre, on travaille le plus souvent en coordonnées relatives, du fait de la difficulté à évaluer quels seront les
points pouvant être considérés comme "fixes", et aussi de la difficulté qu'il y aurait à déterminer pour eux des coordonnées de départ aussi précises que pourront être localement les mesures. Dans la pratique on prendra donc le plus souvent un seul point ( ou éventuellement un sous ensemble de points) dont on fixera arbitrairement les coordonnées de départ, et on fixera une direction arbitraire, éléments par rapport auxquels ont calculera l'ensemble des autres coordonnées relativement à celles ci.

Dans le programme "TRIANGU" on calcule des azimuts et des distances entre points géodésiques connus en coordonnées (voir **Problème Direct**) ce qui revient à des calculs de géodésiques

Le rayon de courbure d'une géodésique est (Formule d'Euler):

$$ (2) \qquad \frac{1}{R} = \frac{\cos^2 Az}{\rho} + \frac{\sin^2 Az}{N} $$

avec $\rho$ rayon de courbure de l'ellipse méridienne

$$ \rho = \frac{a \cdot (1 - e^2)}{(1 - e^2 \sin^2 \Phi_M)^{\frac{3}{2}}} \quad \text{et M point "moyen" de la géodésique} $$

donc

$$ \frac{1}{R} = \cos^2 Az \left[ \frac{(1 - e^2 \sin^2 \Phi_M)^{\frac{3}{2}}}{A \cdot (1 - e^2)} \frac{\sqrt{1 - e^2 \sin^2 \Phi_M}}{a} + \frac{\sqrt{1 - e^2 \sin^2 \Phi_M}}{a} \right] $$

soit

$$ \frac{1}{R} = \frac{1 + e^2 \cos^2 Az \cos^2 \Phi_M \sqrt{1 - e^2 \sin^2 \Phi_M}}{a} $$

donc (3)

$$ \frac{c}{\sqrt{1+e^2 \cos^2 \Phi_M} \cdot [1+e^2 \cos^2 \Phi_M \cos^2 Az]} $$

avec

$$ c = \frac{a^2}{b} \quad \text{et} \quad e^2 = \frac{e^2}{1 - e^2} $$
Réduction des mesures à l'ellipsoïde :

Pour simplifier les calculs, on assimile localement la géodésique à un arc de cercle dont le rayon est le rayon de la géodésique ce qui permet d'utiliser les formules de trigonométrie sphérique. Cette simplification n'entraine pas d'erreur notable pour des distances entre points allant jusqu'à 60 km.

Le passage de la géodésique à l'arc de cercle permet d'utiliser les formules habituelles de trigonométrie sphérique pour passer des éléments angulaires aux éléments métriques.

On donne ci-dessous les principales formules utilisées pour réduire les données à l'ellipsoïde :

Distances:

<figure><img src="image_56.png" /><figcaption>Figure.2</figcaption></figure>

Dans le triangle sphérique PA $oB_0$, on calcule la longueur de l'arc de sphère $A_0B_0$ (géodésique) par (fig.2):
$$ (4) \qquad \cos A_0 B_0 = \sin\Phi_B \cdot \sin\Phi_A + \cos\Phi_B \cdot \cos\Phi_A \cdot \cos(\lambda_B - \lambda_A) $$

soit:

$$ (5) \qquad A_0 B_0 = R \cos^{-1} \left[ \frac{\sin Y_A \cdot \sin Y_B}{R^2} + \frac{\cos Y_A \cdot \cos Y_B \cdot \cos(X_B - X_A)}{R^2} \right] $$

### Angles horizontaux (Fig.3):

Pour simplifier la formulation, dans la figure 2, on appelle A, B, C les projections des points d'observation correspondant sur la sphère de rayon R. Dans le triangle sphérique PAB,

$$ (6) \qquad \cos PB = \cos AB \cdot \cos PA + \sin AB \cdot \sin PA \cdot \cos G_{AB} $$

$$ (7) \qquad \cos G_{AB} = \frac{\cos\left[\frac{\pi}{2} - \Phi_B\right] \cdot \cos AB \cdot \cos\left[\frac{\pi}{2} - \Phi_A\right]}{\sin AB \cdot \sin\left[\frac{\pi}{2} - \Phi_A\right]} = \sin\frac{Y_B}{R} \cdot \frac{\cos\frac{AB}{R} \cdot \sin Y_A}{\sin\frac{AB}{R} \cdot \cos Y_A} $$

et (8)

$$ \sin G_{AB} = \frac{\sin PB \cdot \sin P}{\sin AB} = \frac{\sin \frac{X_B - X_A}{R} \cdot \cos \frac{Y_B}{R}}{\sin \frac{AB}{R}} $$

il vient (9)

$$ \tan G_{AB} = \frac{\sin \frac{X_B - X_A}{R} \cdot \cos \frac{Y_B}{R} \cdot \cos \frac{Y_A}{R}}{\sin \frac{Y_B}{R} \cdot \cos \frac{AB}{R} \cdot \sin \frac{Y_A}{R}} \quad t $$

permettant de calculer le gisement G<sub>AB</sub>. De la même façon on calcule le gisement G<sub>AC</sub>, et l'on obtient l'angle horizontal au point A:

$$ (10) \qquad \alpha = C_{AB} = G_{AB} - G_{AC} $$
<figure><img src="image_57.png" /><figcaption>Figure 3.</figcaption></figure>

### Passage des éléments géométriques mesurés aux éléments sphériques correspondants:

*Passage de la distance entre repères à la corde:

$$ (11) \qquad D_0 = \sqrt{\frac{(D_{\pi}^2 - D_z^2) \cdot R}{(R+Z_A)(R+Z_B)}} $$

*Passage de la corde à l'arc:

$$ (12) \qquad D_{\text{arc}} = 2 \cdot R \cdot \arctan \left[ \frac{D_0}{\sqrt{4R^2 - D_0^2}} \right] $$

R étant le rayon de courbure de la géodésique.

## LA PROJECTION UTM ET LA REDUCTION DES MESURES A CETTE PROJECTION

Le système de projection U.T.M. (Universal Transverse Mercator) est un système de projection conforme consistant en un assemblage de fuseaux occupant chacun 3
degrés de longitude par rapport à un méridien origine choisi sur une zone déterminée. Chaque fuseau couvre l'ensemble des latitudes (d'un pôle à l'autre). L'image du méridien origine est une droite parallèle à l'axe des ordonnées. L'image de l'équateur est une droite perpendiculaire. Ce système permet de remplacer les calculs complexes sur l'ellipsoïde par des calculs très simples dans le plan moyennant deux types de corrections :

- une correction d'échelle E, définie par:

$$E = \frac{ds}{DS}$$

qui est le rapport de la longueur de l'élément d'arc plan à l'élément d'arc ellipsoïdique homologue (le facteur d'échelle est égal à 0.9996 sur le méridien origine, et variable selon la position de l'arc en dehors). - une correction angulaire de réduction à la corde : l'image ab d'une géodésique AB sur l'ellipsoïde n'est pas le segment ab, mais un arc de courbe ab. En outre on doit faire intervenir la convergence des méridiens c de la projection; l'azimut Z d'un segment sur l'ellipsoïde est relié au gisement G du segment homologue sur la projection par:

$$Z = G + c - p_a, \text{ avec}$$

Z = angle entre le Nord géographique et la direction AB,

G = angle entre l'axe des ordonnées de la projection au point a et la direction du segment ab,

c = angle entre la direction du méridien passant par a dans la projection et celle de l'axe des ordonnées,

p_a= angle entre la direction du segment ab sur la projection et celle de l'arc ab (image de AB), au point a.

Des formules relativement simples d'emploi permettent d'effectuer ces corrections et sont programmées comme des "fonctions" dans TRIANGU.

Le passage, direct et inverse, des coordonnées géographiques aux coordonnées UTM, s'effectue par des sous-programmes usuels appelés respectivement GEOUTM et UTMGEO.

## PROBLEME DIRECT ET DERIVEES PARTIELLES

Les éléments géométriques utilisés sont les distances entre repères et les angles horizontaux.
Distances:

A partir des coordonnées X, Y, Z des points A et B, on calcule successivement les éléments A₀, B₀, D₀ et AB.

Calcul des dérivées partielles de AB:

Dans le programme TRIANGU les altitudes des points sont considérées comme parfaitement connues et par conséquent invariables. On a donc :

$$ (13) \qquad \frac{\partial AB}{\partial X_B} = \frac{X_B - X_A}{AB} \qquad \frac{\partial AB}{\partial X_B} = \frac{X_B - X_A}{AB} $$

$$ \frac{\partial AB}{\partial Y_A} = \frac{Y_B - Y_A}{AB} $$

$$ \frac{\partial AB}{\partial Y_B} = \frac{Y_B - Y_A}{AB} $$

Angles horizontaux:

Le problème direct est exprimé à partir des coordonnées des points A, B, C par les formules (9) et (10)

$$ (14) \qquad \alpha = C_{AB} = G_{AB} - G_{AC} $$

Pour le calcul des dérivées partielles, on considérera que l'on peut assimiler l'angle dièdre en A à celui du cas plan; $\alpha_0 = \tan^{-1} G_{AB} - \tan^{-1} G_{AC}$. Alors les dérivées partielles de $\alpha$ s'écrivent:

$$ (15) \qquad \frac{\partial \alpha}{\partial X_A} = \frac{Y_B - Y_A}{AB_0^2} + \frac{Y_C - Y_A}{AC_0^2} $$

$$ \frac{\partial \alpha}{\partial X_B} = \frac{Y_B - Y_A}{AB_0^2} $$

$$ \frac{\partial \alpha}{\partial X_C} = \frac{Y_C - Y_A}{AC_0^2} $$

$$ \frac{\partial \alpha}{\partial Y_A} = \frac{X_B - X_A}{AB_0^2} + \frac{X_C - X_A}{AC_0^2} $$

$$ \frac{\partial \alpha}{\partial Y_B} = \frac{X_B - X_A}{AB_0^2} $$
$$ \frac{\partial \alpha}{\partial Y_C} = \frac{X_C - X_A}{A C_0^2} $$

## MISE EN OEUVRE DE LA PROCEDURE D'INVERSION : LE PROGRAMME TRIANGU

Le logiciel TRIANGU se compose d'un programme principal (Fig.4) traduisant la mise en oeuvre du processus d'inversion (Fig.5) et de 19 sous-programmes ou fonctions pour le traitement des entrées-sorties et d'algorithmes particuliers du processus d'inversion.

### entrées :

Trois fichiers permettent d'introduire les différents types d'information nécessaires au traitement

-Les paramètres de la projection (.PRJ) : Le rayon polaire de l'ellipsoïde (c), la seconde excentricité, le numéro de fuseau de la projection UTM.

-Les coordonnées à priori des stations et points (.COR) et leurs domaines de variation. Cette entrée est gérée par le sous-programme INPUT1.

-Les données mesurées et leur intervalle de confiance par le sous-programme ENTDON. La nature de la donnée est identifiée par un code IDON permettant de sélectionner le type de traitement adapté (sur les angles horizontaux ou sur les distances).

### sorties :

Elles sont gérées par le programme principal pour les sorties auxiliaires de contrôle ou par le sous-programme RESULT pour les sorties finales sur imprimante ou fichier.

Une option permet une sortie des coordonnées compensées finales dans les différents systèmes de coordonnées : UTM, géographiques (en grades, degrés décimaux ou degrés minutes, secondes) ou cartésiennes géocentriques.

### Module de calcul du problème direct :

Ce sous-programme PBDIR a pour objet de déterminer les valeurs calculées à partir des valeurs courantes des paramètres pour chaque observation ainsi que les dérivées
partielles correspondantes. Des tests (sur IDON) permettent de sélectionner les traitements spécifiques pour chaque type de données.

### Modules relatifs à l'inversion :

Deux modules permettent de calculer :

a) le second membre de la relation de récurrence

$$P_{k+1} = P_k + (G_k^T C_d^{-1} G_k + C_p^{-1})^{-1} (G_k^T C_d^{-1} \{d_0 - g(p_k)\} - C_p^{-1} (p_k - p_0))$$

par le sous programme FSGM

$$S = G_k^T C_d^{-1} [d_0 - g(p_k)] - C_p^{-1} (p_k - p_0)$$

b) la matrice de résolution par le sous-programme MATRES

$$A = G_k^T C_d^{-1} G_k + C_p^{-1}$$

On résoud ensuite le système linéaire $A \cdot X = S$ par le sous-programme DCHOL qui met en oeuvre cette résolution par la méthode de Choleski, et fourni un vecteur X dont les éléments sont les corrections à apporter aux paramètres à l'issue de la kième itération.

Le test de convergence du processus itératif est effectué par la fonction SIGMA dont l'expression est de la forme :

$$\sigma = \sqrt{\frac{1}{2n} \sum \left[ \frac{d_0 - d_c}{E_d} \right]^2 + \frac{1}{2m} \sum \left[ \frac{p_0 - p_k}{E_p} \right]^2}$$

avec $d_c = g(p_k)$ valeur calculée des données à partir des valeurs à l'itération k des paramètres; $E_d$ l'intervalle de confiance de la donnée; et $E_p$ l'écart a priori sur le paramètre; les sommes étant faites sur les n données et les m paramètres.

Les résidus sur chaque donnée sont sortis sur option ou après arrêt du processus itératif; lorsqu'il s'agit d'une donnée angulaire on calcule également l'équivalent en distance de ce résidu :
$r_d = d'_o - d'_c = (d_o - d_c) \cdot \Phi_{AB}$ c'est-à-dire la déviation métrique de l'écart angulaire correspondant à la distance sur lequel il s'applique.

On calcule par ailleurs, à chaque itération l'erreur quadratique moyenne des résidus "ramenés aux distances" par la fonction RESID:

$$ \sigma_{rd} = \sqrt{\frac{\sum_{n} \left[ \frac{d'_{o} - d'_{c}}{E_{d}} \right]^{2}}{n}} $$

Chaque donnée est traitée différemment selon qu'il s'agit d'une donnée métrique (distance ou dénivellée) ou angulaire (angle azimutal ou zénithal) :

- donnée métrique : $d'_{o} - d'_{c} = d_{o} - d_{c}$

- donnée angulaire : $d'_{o} - d'_{c} = (d_{o} - d_{c}) \cdot \text{dist.}$

En cas de divergence (problème non linéaire mal contraint), il se peut qu'à la kième itération $\sigma_k > \sigma_{k-1}$. Dans ce cas on divise par 2 le facteur correctif X puis on repart dans le processus itératif avec les nouveaux paramètres ainsi calculés, puis on teste à nouveau sur $\sigma$, et ainsi de suite jusqu'à obtenir la convergence. Ce processus est caractérisé par la valeur du paramètre ALPHA dans le programme principal.

Le calcul s'arrête soit en fonction du nombre d'itérations demandé par l'opérateur, soit lorsque $\sigma_{rd}$ devient inférieur à une valeur arbitraire fixée à l'avance. A la fin du processus itératif on calcule la matrice de variance-covariance a posteriori des paramètres en inversant la matrice de résolution.

$$ M_k = G_k^T C_d^{-1} G_k + C_p^{-1} $$

Cette matrice est inversée par la méthode de Choleski par le sous-programme DCHOLI. Les erreurs a posteriori sur les paramères finaux sont alors obtenu par la racine carrée des éléments diagonaux de $M^{-1}$.

Le calcul des valeurs propres de $M^{-1}$ permet alors d'obtenir les valeurs principales et les orientations des axes des ellipsoides d'erreur sur chaque point.

Les autres sous-programmes ou fonctions concernent des calculs proprement géodésiques.
<figure><img src="image_58.png" /><figcaption>Figure.4 - Schéma de principe du programme TRIANGU</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure.5 - Algorithme général d'inversion d'après la méthode Tarantola & Valette</figcaption></figure>
BIBLIOGRAPHIE

- Tarantola A. et B. Valette, Inverse Problems = Quest for information, J.Geophys., 50,pp.159-170, 1982.

- Tarantola A. et B. Valette, Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. Space Phys., 20,no 2, pp.219-232, 1982.

- Dupuy M. et Dufour H.M., La géodésie, Collection Que sais-je, Presse Universitaires de France, 1969.