Brevet FR2778474A1 — publié le 19991112

-1- 2778474 "Procédé pour reconstruire une surface avec mise en oeuvre d'une prédiction linéaire et applications de ce procédé" La présente invention concerne un procédé pour reconstruire une surface avec mise en oeuvre d'une prédiction linéaire. Elle vise également des applications de ce procédé. Cette invention s'inscrit dans le cadre du développement et de l'application de méthodes de statistique et de traitement du signal aux techniques de la mesure. Elle est basée sur la notion classique de prédiction linéaire d'une fonction aléatoire Z, qui consiste à extraire le maximum d'information des mesures effectuées en des sites S, et en particulier d'estimer une densité de probabilité du processus mesuré, en des points X situé en dehors des points mesurés. La prédiction linéaire d'une fonction aléatoire Z, consiste à calculer sa densité de probabilité conditionnelle, conditionnée par les valeurs de Z connues par mesurage en un certain nombre de points. L'intérêt est d'estimer ses valeurs probables en tout point o Z n'a pas été mesurée. La prédiction linéaire est utilisée dans de nombreux domaines des sciences: biologie, économétrie, sciences sociales, physique appliquée, géostatistique... La prédiction temporelle concerne celle de fonctions du temps, en général en extrapolation (cours en bourse, évolution de la démographie etc.). La prédiction spatiale vise l'affinage de l'information apportée par un échantillonnage d'une fonction dont la variable est spatiale (dimension 1 à 3). Ces aspects sont notamment décrits dans les références bibliographiques [1], [4], [9], [11]. En ce qui concerne la métrologie des surfaces fabriquées en génie mécanique et génie optique, la mise -2- 2778474 en oeuvre de la prédiction a fait l'objet de travaux récents effectués par le présent déposant ayant fait l'objet des publications [5] à [7]. Les figures suivantes sont représentatives de l'état de la technique dans ce domaine: - la figure 1 représente un défaut de forme (en mécanique) ou "Pic à Vallée" (en optique); - la figure 2 illustre un exemple d'évolution de l'espérance conditionnelle Ec(X); - la figure 3 illustre un exemple d'évolution de la variance conditionnelle Vc(X); - la figure 4 illustre un exemple d'évolution des bornes des intervalles de confiance de (Z(X) IZ(S)); - la figure 5 représente une variation de la prédiction d'un profil 2D; - la figure 6 est une représentation graphique d'un fichier de 10OxlOO points mesurés; - la figure 7 est une représentation graphique d'un fichier échantillonné 33x33 puis prédit. La métrologie des surfaces fabriquées industrielles est surtout orientée vers le contrôle des spécifications. Celles-ci sont en général constituées de limites imposées à certaines grandeurs de la surface. Pour la plupart, ces grandeurs seraient nommées, en mathématiques, des normes '1(1 1, II 112)' En mécanique, la norme I II est le défaut de forme au sens de la norme ISO 1101, en référence à la figure 1. En rugosimétrie des pièces mécaniques, ce défaut est le paramètre d'indice t: Pt, Rt. En optique, il est appelé " défaut Pic à Vallée " (P-V) ou plus souvent Peak to Valley (PTV). La norme H H2 est définie comme l'écart quadratique moyen par rapport à une surface moyenne, en général la surface des moindres carrés estimée sur les points 3- 2778474 mesurés. Appelée ici " modèle ", elle dépend d'un paramètre vectoriel p, et notée Mp. Cet écart quadratique moyen est le "Mean Square", sa racine carrée est le "Root Mean Square" (RMS) c'est à dire l'écart-type. C'est une spécification usuelle pour les surfaces; en rugosimétrie des surfaces mécaniques, il fournit les paramètres d'indice q: Pq, Rq. On connaît déjà une méthode d'estimation du défaut de forme à partir des bornes des intervalles de confiance de la surface, ce qui ne permet pas de déterminer les intervalles de confiance du défaut de forme. La référence bibliographique [7] propose une estimation du modèle des moindres carrés Mp à partir de toute la densité de probabilité conditionnelle, c'est à dire de (ZxlZs) et pas seulement partir de Zs comme le font toutes les méthodes classiques. On peut ensuite estimer l'espérance et la variance du MS conditionnées par Zs. La précision obtenue sur l'estimation du MS est grandement améliorée par rapport aux méthodes non statistiques qui ne peuvent prendre en compte que les données Zs, et ne fournissent pas de variance (c'est à dire pas d'incertitude). Soit Z une fonction aléatoire gaussienne stationnaire, de moyenne nulle. Z est à valeurs dans Rp. Z dépend d'une variable "événement" o appartenant à un ensemble Q de possibles, ainsi que d'une variable "spatiale" x appartenant à un domaine Dc R". Par exemple, Z représente une surface d'une pièce mécanique fabriquée en série. o est le numéro de la pièce dans la série. La variable spatiale est un point (Xy) de R2 et D représente les limites matérielles de la surface. Z(Xy) est l'altitude de la surface au point considéré. -4- 2778474 Z est mesurée en S points de D, appelés ici sites de mesure, et S=(S1,...,Ss) l'ensemble de ces sites. Soit X=(X1,---,Xx) un ensemble de X points o Z est inconnue et doit être estimée. Pour alléger les notations, on écrira5 "Zx" pour Z(x). Appelons Cor:R ±>R, sa fonction d'auto-corrélation. On rappelle que la covariance de deux vecteurs aléatoires réels A et B est Cov(A, B)=E[A.tB], o E désigne l'espérance et t la transposée. Czzx =Cov(Zx Zx)=(E[Zx i.Zx i est la covariance (non conditionnelle) de Z sur l'ensemble X. On définit les vecteurs et matrices de covariance suivants: Cz-zs= Cov(Zx Zs)=(E[ZX [zi]) _A CzZ.x=cov(zszx) =(E[Zs =.zX J]).- = tczxzs CZszs =Cov(zs Zj)=(E[ ZSi zs ji -iSj lS(équations (1.1) Un résultat classique de mathématique, illustré par la référence [ 3] pp. 69 sq.) énonce que les deux premiers moments conditionnels du vecteur aléatoire inconnu Zx(z8, x) conditionné par le vecteur X) aléatoire connu Zs=(zs,..,Zs), sont l'espérance conditionnelle Ec(X): Ec(X)=E[ZxlZS]=C.C -.Zs (1.2) zxzs zs,zs et la variance conditionnelle Vc(X): Vo(x)=V[ZxIZs]=E[lIzs]=Czxzx-c C6.C (1.3) zx,zs zsszs zs,Zx Le vecteur aléatoire Zx étant gaussien, sa densité de probabilité est entièrement déterminée par ces deux premiers moments. Dans la suite, " prédire Zx " aura le sens " exprimer ces deux moments conditionnels " (1. 2) et (1.3), définis par les relations de covariance (1.1). Il est à noter que les équations (1.2) et (1.3) ne sont pas définies sans les relations de covariance (1.1), qui en constituent le coeur. La figure 2 et la figure 3 illustrent le résultat de la prédiction d'une fonction réelle de deux variables réelles connue par ses valeurs aux points d'un maillage carré (0,1,2,3) et inconnue ailleurs. L'espérance conditionnelle est une interpolation qui passe effectivement par les points mesurés. Cette espérance Ec(X) représente la position la plus probable de Z(X). En cela, c'est une " reconstruction " de la surface, meilleure qu'une interpolation polynomiale arbitraire. La variance conditionnelle Vc(X) exprime l'incertitude sur la position de Z(X) par rapport à Ec(X)- A partir de ces deux premiers moments conditionnels qui définissent entièrement la variable aléatoire gaussienne Z(X), on peut déduire des intervalles de confiance Ec(X) ka. (X), qui contiennent Z(X) au taux de confiance a, en référence à la figure 4. Comme premier exemple, la figure 5 représente la prédiction d'un relevé rugosimétrique 2D d'une pièce mécanique, connue par un échantillonnage de 9 points. Comme second exemple, voici un fichier de 10OxlOO0 points mesurés sur un miroir avec un interféromètre Zeiss -6- 2778474 à contraste de phase (Figure 6). On en a extrait 33 x 33 points, seuls connus de la méthode prédictive. La figure 7 est la reconstruction de la surface par prédiction à partir de cet échantillonnage du fichier de référence. Le but de la présente invention est de proposer un procédé de reconstruction plus performant que les procédés d'estimation actuels. Suivant l'invention, le procédé pour reconstruire une surface connue par un ensemble de mesures en des sites, mettant en oeuvre une prédiction linéaire en des points situés en dehors des sites, est caractérisé en ce qu'il comprend une utilisation de mesures de variables géométriques associées à ladite surface, et en ce qu'il procure une estimation optimale de caractéristiques desdites variables géométriques. Il peut s'agir soit de mesures d'altitudes, soit de mesure de pentes, avec dans le premier cas, une estimation optimale des caractéristiques géométriques de la surface, et dans le second cas, une estimation des caractéristiques des pentes. - Lorsque ce procédé est mis en euvre pour la prédiction de Z(X) conditionné par des altitudes Z(S) mesurées, il comprend en outre un affinage de l'estimation de la fonction d'auto- corrélation et de la fonction de variogramme de la fonction aléatoire Z, cet affinage comprenant une utilisation de la prédiction effectuée à partir d'une première estimation grossière desdites fonctions d'auto-corrélation et de variogramme. Le procédé de reconstruction selon l'invention peut avantageusement être mis en ouvre pour la prédiction des altitudes Z(X) conditionnées par des pentes Z' (S) mesurées. Dans cette version de l'invention, on adapte aux pentes, selon la référence bibliographique [7], une solution des moindres carrés généralisée: 7. 2778474 p = (tfQ.fy.) t f.Q. y Dans une mise en oeuvre du procédé selon l'invention pour la prédiction des pentes Z'(X) conditionnées par des pentes Z'(S) mesurées, celui-ci comprend en outre une détermination d'un vecteur complet des pentes mesurées et prédites. Une première application du procédé de reconstruction selon l'invention concerne la rugosimétrie tridimensionnelle de pièces mécaniques. Dans cette application, on réalise une reconstitution visuelle par tirage aléatoire d'une réalisation de Z(X) avec comme densité la densité de probabilité conditionnelle de (ZxlZs). Un second domaine d'application du procédé de reconstruction selon l'invention concerne la métrologie tridimensionnelle de pièces mécaniques. Une autre voie d'application préférée du procédé selon l'invention concerne sa mise en ouvre dans un système d'analyse de surface d'onde par interférométrie, ce système comprenant une analyse de réseaux de franges d'interférences conduisant à -un nuage de points tridimensionnel S, un traitement de ces points S par prédiction pour la reconstruction d'un front d'onde et une estimation des caractéristiques de ce front d'onde. Le traitement mis en oeuvre dans ce procédé comprend alors: - une égalisation d'histogramme adaptative, - Une squelettisation des franges, une recherche des niveaux, - une recherche optimisée des extrema des franges, et - un traitement par prédiction. -8- 2778474 D'autres particularités et avantages de l'invention apparaîtront encore dans la description ci-après. Aux dessins annexés donnés à titre d'exemples non limitatifs: - la figure 8 représente un maxi et un mini de l'espérance conditionnelle; - la figure 9 représente un exemple de disposition des points de simulation; - la figure 10 illustre le concept de pente maximale sur l'espérance conditionnelle; la figure 11 illustre une disposition des points de simulation; - la figure 12 illustre une prédiction de Z(X) conditionné par Z'(S); - la figure 13 illustre une prédiction de Z(X) conditionné par Z'(S); - la figure 14 illustre une disposition des points pour la simulation d'un pavé n 1; - la figure 15 illustre une disposition des points pour la simulation d'un pavé n 2; - la figure 15 illustre un exemple de relevé rugosimétrique tridimensionnel 5x10 mm: 100 profils espacés de 0,05 mm; - la figure 17 illustre un exemple de relevé rapide, avec 17 profils espacés de 0,3 mm; la figure 18 illustre un exemple de reconstruction par prédiction, avec 100 profils; - la figure 19 représente un cercle usiné, mesuré sur machine de circularité de type "Talyrond"; - la figure 20 illustre une comparaison de l'ajustement d'un modèle des moindres carrés avec et sans prédiction; la figure 21 illustre une comparaison des estimations avec et sans prédiction de l'écart quadratique moyen (MS); 9 2778474 - la figure 22 illustre une comparaison des estimations avec et sans prédiction du défaut de forme macro-géométrique (PTV); - la figure 23 représente un exemple de réseau complexe de franges d'interférences; - la figure 24 représente un exemple de squelettisation des franges sombres; la figure 25 illustre une technique d'insémination; - la figure 26 illustre un exemple d'ajustement d'un polynôme à une frange claire pour la détermination du sommet; - la figure 27 illustre un exemple de disposition des points de prédiction X pour la commande par la méthode de Shack-Hartmann; et - la figure 28 illustre une comparaison des erreurs de reconstruction montrant les améliorations procurées par le procédé de reconstruction selon l'invention. On va maintenant décrire plusieurs exemple de mise en oeuvre du procédé de reconstruction selon l'invention, en référence aux figures 8 à 28. On se propose tout d'abord de mettre en oeuvre une prédiction de Z(X) conditionné par les altitudes Z(S) mesurées. On propose ici d'utiliser une méthode expérimentale d'estimation de la densité de probabilité de grandeurs (telles que PTV ou RMS) issues d'une réalisation de la fonction aléatoire Z. C'est une méthode de type Monte- Carlo, consistant en un grand nombre de tirages au hasard O de " trajectoires " Z(o). En fait, on se limite à un tirage de vecteurs aléatoires Z(X)=(Z(X1),-.,Z(X)), o X est le plus grand nombre possible, compatible avec les moyens de calcul (capacité mémoire et temps de calcul). Dans les relations (1.2) et (1.3), les matrices d'indice -10- 2778474 X ont un nombre de lignes ou de colonnes égal à X. Si X = 1000, Czxzx est une matrice 1000xl00. Or 1000 points forment un quadrillage 33x33, ce qui est peu. D'o la nécessité d'une répartition judicieuse des points. Le5 tirage aléatoire de Z(X) se fait en tirant une réalisation d'une variable aléatoire vectorielle gaussienne de dimensionX, de matrice de covariance Cz$,zs(que l'on ne calcule qu'une seule fois), et de moyenne Eó(X) En référence à [7], on définit g comme le signal résiduel g(X)=Z(X)-E (X). Sur la figure 8, on a représenté les zones au voisinage du maxi et du mini de l'espérance conditionnelle. C'est dans ces zones que la probabilité d'extrema de Z est la plus forte, et que l'on choisit les points X de simulation de Monte-Carlo, comme indiqué sur la figure 9, et on répartit les points de simulation X=(X1,.-.,Xx) simultanément sur les deux zones 1 et 2. Pour chaque tirage aléatoire (Z(X1),--.,Z(XX))IZs, on obtient un PTV du tirage: PTV= max Z(Xi)- min Z(Xi). La série statistique de ces PTV donne un histogramme H. Une variante consiste à simuler séparément en zone 1 (donnant un histogramme H1) et en zone 2 (donnant un histogramme H2). Comme PTV=max-min, l'histogramme H du PTV est la convolution H=H1*H2. Ceci suppose que les zones 1 et 2 soient assez éloignées pour que Z dans l'une et dans l'autre soient décorrélés. On en tire un histogramme représentant la densité de probabilité de la variable aléatoire PTV. Cette densité fournit une espérance (la valeur estimée) et des intervalles de confiance. En référence à [7], la décomposition de Z(X)=Ec(X)+g(X) permet d'écrire: -11- 2778474 MSw= f(Z2(x)lZ(S))dX = (Ec (X) + g( X d xC 2e = (Ec (X))2.dX + 2. Er (X).g(X)d+ (g (X)).d2 Comme le signal résiduel g est " tiré au hasard " dans la simulation, c'est lui qui représente la composante aléatoire de la fonction a posteriori (ZxIZs). Pour chaque (, on obtient une " trajectoire " go et donc une valeur MS. On en tire un histogramme représentant la densité de probabilité de la variable aléatoire MS. Ici également, cette densité fournit une espérance (la valeur estimée) et des intervalles de confiance. Comme la10 densité est connue, on peut aussi en déduire des intervalles de confiance du RMS qui est la racine carrée du MS. La difficulté réside dans la non-stationnarité du terme Ec (X).gw (X), imposant de simuler X sur tout le domaine D, avec un gros volume de calcul. Si on veut calculer la variance du MS, non à partir d'une simulation de Monte-Carlo,, mais selon la référence [7], on se heurte à la méconnaissance de la covariance de g2 Des approximations, bien confirmées par expérimentation sur des miroirs, conduisent à: 2. v VarMS=4 Ec.Vc o E 2 Ec2z(u).dX et Vc=u Vcz(u).dX sont les c = fueD fu eD moyennes respectivement de la variance et du carré de l'espérance conditionnelle. On va maintenant, dans le cadre du procédé de reconstruction selon l'invention, estimer la fonction d'auto-corrélation et le variogramme. On rappelle les définitions classiques: -12- 2778474 Si Z est stationnaire, sa covariance Covz(X,X+T)=E[Zx. Zx+T] est indépendante de X et ne dépend que du " pas " T; elle prend alors le nom de fonction d'auto-corrélation et est notée ici Corz(T). Par ailleurs, le variogramme Varz(T)=E[(ZX+T -Zx)] est la fonction de variance des incréments au pas T (ou " fonction de structure "), et est reliée à la fonction d'auto-corrélation Cor par la relation: Varz (T) =2 (Corz (0) - Corz (T)). Les notions de fonction d'auto-corrélation (probabiliste) et de d'auto-corrélation spatiale se confondent sous l'hypothèse d'ergodicité, qui est expérimentalement bien vérifiée si Z est stationnaire, et en pratique si les basses fréquences sont retirées avant prédiction. Par la suite, on définit donc Corz(T) comme la moyenne spatiale: xCorz(T)={ Zx.ZX+Tdx JxD et Varz(T) comme: Varz) (ZXT ZX)dx Var z(T)= ZX+TZ)dl EXD Corz(T)et Varz(T) sont donc des variables aléatoires qui dépendent du " tirage " O et dont on cherche à exprimer les densités de probabilité, ou au moins les deux premiers moments. Par la décomposition Z(X)=Ec(X)+g(X), l'espérance conditionnelle de la variable aléatoire Corz(T) est: cCorz(T)= Ecf ZX. ZX+TdX =EC f (Ec (X)+g (X)).(Ec (X +T) + g (X + T)) d XED xED -13- 2778474 Après calculs: cCor z (T) = Ec (X).Ec (X + T) + g(X).g(X + T) c'est à dire la somme des fonction d'auto-corrélation spatiales de Ec et de g. EcCor z (T) = Cor Ec (T) + Cor g (T) Comme g n'est pas stationnaire, et qu'elle n'est de plus pas toujours connue (sauf les cas o l'on peut faire une mesure zoom sur la surface), on peut prendre comme fonction Corg(T) une fonction affine décroissante, égale à Vc(X) pour T=0 et nulle pour T=T0(X), la distance inter-sites locale. En effet, g apparaît comme décorrélée pour T>To (X). La variance conditionnelle de Cor z(T) est: " 2' Vc Cor z (T) = Vc Zx. ZX+T =EcZx Zx+T) - (EC(X. ZX+T)) Les calculs conduisent, après approximation, à majorer cette variance par Var MS qui est aussi Var(Corz(O)). On peut utilement se référer à [7] pour le calcul de Var MS. Or ce dernier a pour valeur approchée VarMS=4.Ec. Vc. Les intervalles de confiance de Corz(T) sont donc définis par: EcCor z (T) = Cor E (T) + Cor g (T) Vc Cor z (T) 4. Ec_ que l'on utilise dans l'inégalité de Bienaymé- Tchebytchev. L'estimation du variogramme Varz (T)=2.(Corz(0)-Corz(T)) s'en déduit immédiatement. En particulier, la variance de Varz(T) est inférieure à 8.E2c.Vc et son espérance est Ec Var z (T)= Var Ei (T) + Var g (T). -14- 2778474 On va maintenant estimer le PTV et le RMS des "pentes au pas T". En optique, l'une des spécifications sensibles est le PTV des pentes au pas T, c'est à dire des incréments (ou différences finies) de Z au pas T: ÈZT (X)=ZX+T-ZX. Le MS de ces pentes n'est autre que Varz(T)/T2 que l'on estime comme précédemment. On choisit les sites S et les points de simulation X dans les zones 1 et 2 définies sur les figures 10 et 11 de telle sorte que la plus grande pente de l'espérance conditionnelle s'y observe. La simulation de Monte-Carlo de ce PTV peut se faire par la méthode suivante: Pour chaque tirage aléatoire o des (ZxIZs), on prend toutes les différences fines IZxm -ZXnI telles que IlXm-Xn |=T, puis le maximum de ces différences. De ces maxima, on tire un histogramme du PTV des pentes qui converge vers la densité de probabilité du PTV. On en déduit l'espérance et des bornes des intervalles de confiance du PTV des pentes du signal résiduel g. On rajoute alors la plus grande. pente de l'espérance conditionnelle. On va maintenant prédire les altitudes Z(X) conditionnées par les pentes Z'(S) mesurées. On suppose ici que la technologie de mesurage fournit Z en chaque site de mesure Si(xy) les pentes Z;= - (y) et Z2=a (y). Z est ici supposée définie sur Dc:R2, à valeurs dans R, et de classe C1. Dans des situations pratiques dans lesquelles Z peut n'être dérivable en aucun point comme une fractale, on considère le phénomène physique réel, filtré passe-bas - 2778474 par le procédé de mesurage, de sorte que les dérivées partielles soient définies. Dans l'hypothèse o l'on s'intéresse à la densité de probabilité du signal non filtré, il faut appliquer la démarche suivante: on effectue une convolution de la densité de probabilité issue de prédiction avec l'histogramme expérimental du signal de haute fréquence, obtenu par un zoom (rugosimétrie, ou une sous-pupille en optique). Il s'agit maintenant de prédire Z(X) à partir des pentes Z' (S) mesurées. Définissons les vecteurs aléatoires supposés gaussiens: Z'(s)= (Z;(S1),Z2(S1),...,Z;(S),Z'(S)) et Z(X)=(Z(X),,Z(X Considérons le vecteur aléatoire: (Z'(S)) et exprimons la densité de probabilité de Z(X) conditionnée par Z'(S). Les relations (1.2) et (1.3) s'appliquent: E[ZxlZs]=C.C.Z4 (1.4) ZX. Zs %s % Zx z',z - et la variance conditionnelle Vc(X): V[ZxIZ]=E[z2 lz]=Czx C czx-1 c. (1-5) Xzx,z's, zs zs,zx Ces relations ne seraient en rien novatrices si elles se suffisaient à elles-mêmes. Mais sans les relations précises de covariance, elles sont vides de sens. Les termes "C" de (1.4) et (1.5) sont définis de la façon décrite ci-après. Une démarche appropriée est de choisir une fonction d'auto-corrélation de Z (notée comme précédemment Corz), et d'en déduire les corrélations entre Z et Z' -16- 2778474 (notée CorzZ,) et entre Z' et Z' (notée Corz,). Au niveau de l'estimation, on cherchera à estimer directement Corzen fonction des données Zs. On sait par la référence [8], tome 2 page 32, que la covariance de la dérivée f'(c,x) d'une fonction aléatoire d'une variable spatiale f(û,x), est la dérivée seconde de la covariance de f: Covf.(u,v) =E(f'(u).f'(v)) = aCvf (u,v) ôuav Pour deux fonctions aléatoires f et g, notons leur covariance Covf.g(Uv)=E(f(u).g(v)). Soient deux points X1=(xl,yl) et X2=(x2,y2) de D. Ici, Z et Z' sont à 2 variables réelles, et les covariances à 4 variables réelles. Par exemple, Covz.z (X, X2) = Covzz ((x1l, y1), (x2, y2)). Après calculs, on trouve les résultats suivants, toutes les covariances étant fonctions des variables ((X1,y1),(X2,y2)): acovz aCovz Covz,z = óx I C vz22z = ay1 Cov COVZ Covz, V 8x2 8Y2 ________ aX Covz2 Covziz =; ov, a2COVx ó2COVz a2COVz CovZ a= 2Covz aY18x2 z aYl-aY2 Équations (1.6) Cor z est la forme prise par Covz dans le cas stationnaire. On remarque que Covz doit être de classe c2 pour que ces dérivées soient définies et continues. Avec ces fonctions de covariance (1.6), les matrices de covariance des équations (1.4) et (1.5) s'expriment par: -17- 2778474 Czx zX = CovZZ (X, X) = (Covz,z (X,,x))i,j=X... Czx, z = COvzz' (X,S) = (Covz,z. (X, Sj))i =1...,, j=l o Covzz,(Xi,Sj) est un bloc (2,1) égal à Cov.z (X, S,) C, z. = Covz,z (S, X) = (Covz., z (Si, Xj)) =...,j=tX o Covz,.z(Si,Xj)est un bloc (1,2) égal à (Covz;.z (Si,Xi) Covz.z (Si,Xi)) Cz,z = Covz,. (S. S) = (Covzz (si, si))j =1_ o Covz7z, (Si Sj)est un bloc (2,2) égal à COVzizi (Si,Sj) C Vz;.Z2 (Si,S,) CCovz zi(SiSj) CoVz,z(Si,Sj1) (Equations (1)) Cette méthode de prédiction reconstitue donc directement Zx à partir des pentes mesurées Zs. Une validation expérimentale pour une fonction Z d'une variable, est aisée à effectuer: On se donne par exemple Z'(S)=(%0,10,-1) pour S = (0,tZ3,4), ainsi qu'une fonction d'auto-corrélation Corz(T) dont l'allure est indiquée en figure 12. La figure 13 résume le test: Z(X) est entièrement estimé (à une constante verticale près) à partir des seules données Z'(S)=(%0,%0,-1) et Corz(T). La courbe Z' (X) est la dérivée de Z(X). Le test consiste à constater que, pour XeS, les valeurs de cette dérivée sont bien égales aux pentes Z'(S) données. Pour Xz(0,1,Z3,4), Z' (X) est une interpolation des pentes Z'(S). Pour l'estimation de Corz (T), la démarche est classique, en référence à [ 2]. On choisit une famille de fonctions d'auto-corrélation paramétrée par un paramètre - 18 - 2778474 vectoriel que l'on estime au maximum de vraisemblance en fonction des données Z'(S). La méthode de Monte-Carlo permet des tirages aléatoires de Z(X) conditionné par Z' (S) à partir des relations (1.4) et (1.5). L'estimation de Cor et Var et de leurs intervalles de confiance est formellement identique à ce qui a été décrit précédemment. En effet, cette estimation ne fait intervenir que la densité de probabilité conditionnelle, dans laquelle (ZxlZs) est remplacé par (ZxlZs). On va maintenant prendre en compte les non- stationnarités. Comme dans la référence [7], on suppose que la fonction étudiée, nommée ici Y, et dont sont issues les pentes mesurées Ys est de la forme Y=Z+Mp o Z est une fonction aléatoire définie comme précédemment, et Mp un modèle déterministe dépendant linéairement d'un paramètre vectoriel p: Mp (X)=f(X).p avec p=t(p1,...,p,) et f (X) =(fi (X),--.f n (X)) - La méthode courante est d'estimer p à partir des pentes mesurées Yà par une " estimation modale ", en référence à [10].. On obtient une solution des moindres carrés: = (t%. f)-l.tf.YS On propose, dans le cadre du procédé de reconstruction selon l'invention, d'adapter aux pentes une solution des moindres carrés généralisée de même forme que [7]: p = (tfv. Q). ltfi Q. Yf Il est à noter que, avec ou sans prédiction, p est linéaire par rapport à Ys, de la forme p=UYà. Alors M (X)=f(X).ó=f(X).U.Y. On pose F(X)=f(X).U. Alors Me (X) est de la forme Mb(X)=F(X).Y o F(X) ne dépend que de f et S, Y9. 2778474 ainsi que de Cor, si on utilise la prédiction comme en référence [7]. On en déduit M (X)=F. Y_ et Mp (S) =F. YS. On définit les résidus de Y' après avoir retranché le modèle: Zs=Y - Mp(s) On prédit ensuite les (ZxlZs), puis on rajoute le modèle pour obtenir: E[Yx IY] =Mp(X)+E[Zx Iz] =Mp(X)+C.CC.-Z1 Zx,Z s z s,z (S =Fx. YS + C.C -1 *(YS - FS. YS) zxz.sZz',z's Zx,Zs z's.z's d'o E[Y YS]=Fx+C c 0 (I-Fs))Y (1.8) zx,zs.z'sZs On va maintenant considérer la prédiction des pentes Z' (X) conditionnées par les pentes Z' (S) mesurées. E[zXIzS]=C.C' -.z (1.9) z - 20- 2778474 o tous les termes Covz,z,(Ai,Bj) sont des blocs (2,2) égaux à Covzz (Ai, Bj) Covz;,z2 (Ai,Bj)1 (1.11) CoVz.zi (Ai,B,) COVz2.z2(AiBi)J On notera les différences par rapport aux équations (1.4), (1.5) et (1. 7). E[YkIY] =M;(X)+E[ZxjtZ] =UM;(X)+C C - Z' z.x,z.s z's,z's =Mp (X) + C.C -Z.(Y - M (s)) zc,zs z, =F'.YS + C.C -1 -(Yà - FS. Y) z'x,z, z's,zs d'o: E[Y:IY]=F +C C1.(I o W est définie par la relation (1.3) et 124 est la matrice identité d'ordre 2.4. On va maintenant présenter plusieurs applications du procédé selon l'invention. Une première application concerne la rugosimétrie tridimensionnelle de pièces mécaniques. La rugosimétrie tridimensionnelle des pièces mécaniques suppose le déplacement mécanique d'un capteur mécanique (à stylet) ou optique (à laser) au dessus de -21- 2778474 chaque point de la surface à mesurer. Le balayage se fait par un relevé de profils 2D parallèles. La durée du mesurage actuelle est de plusieurs dizaines de minutes à plusieurs heures. Ceci interdit le contrôle " en ligne ",5 c'est à dire de chaque pièce à la sortie de la chaîne de production. Par contre la rugosimétrie 3D a des avantages importants par rapport à la rugosimétrie 2D. La prédiction linéaire permet donc de diminuer le temps de mesurage et d'optimiser l'information apportée par la mesure en espaçant les profils mesurés et en reconstituant les profils manquants. On améliore donc: - la visualisation des états de surface l'estimation des paramètres d'Etat de Surface: PtRtPqRq etc. La reconstitution visuelle se fait par tirage aléatoire d'une réalisation de Z(X) avec comme densité, la densité de probabilité conditionnelle de (ZxlZs). Le tirage se fait de proche en proche: On commence par simuler avec des points X disposés à l'intérieur d'un pavé bordé par deux profils mesurés, en référence à la figure 14. Puis on simule dans un pavé n 2, limitrophe du n l, avec comme points-conditions ceux des profils mesurés qui bordent le pavé 2, ainsi que ceux des points de 1, simulés à l'étape précédente, et qui bordent aussi le pavé 2, en référence à la figure 15. La reconstruction visuelle par prédiction est illustrée par la figure 16 représentant un relevé rugosimétrique tridimensionnel de 5x10 mm avec 100 profils espacés de 0.05 mm. La figure 17 illustre un relevé rapide sur la même surface de 17 profils espacés de 0.3 mm, tandis que la figure 18 représente un exemple de réalisation d'une reconstruction visuelle par prédiction. Des validations de l'estimation des paramètres d'état de surface ont été effectuées en prenant un relevé de -22- 2778474 grand nombre de points (par exemple 2562), extrayant un profil sur 6 (temps de mesurage: 10 mn au lieu d'une heure) et en comparant les intervalles de confiance des paramètres avec la valeur de référence du relevé complet. Les résultats sont satisfaisants puisque toutes les valeurs de référence sont dans les intervalles de confiance, au taux de confiance choisi. Une seconde application industrielle du procédé selon l'invention est la métrologie tridimensionnelle de pièces mécaniques. L'ajustement d'un modèle des moindres carrés à un ensemble de points mesuré est beaucoup plus efficace avec prédiction linéaire que sans prédiction linéaire. Par exemple, la figure 20 est un relevé de 1000 points sur un cercle usiné, avec une amplitude des défauts radiaux de l'ordre du pm. On en a extrait n points répartis régulièrement, qui constituent les sites de mesure S fournis à la méthode prédictive. Le cercle des moindres carrés ajusté avec prédiction linéaire selon les résultats de la référence [7], est beaucoup plus proche du cercle "réel" ajusté aux 1000 points, que le cercle des moindres carrés sans prédiction, en référence à la figure 21. Le défaut de forme est également bien mieux estimé avec le procédé de reconstruction selon l'invention, en référence à la figure 22. Le procédé selon l'invention peut également être mis en oeuvre dans des systèmes d'analyse de surface d'onde par interférométrie. On décrit ci-après une analyse de réseaux de franges d'interférences simples ou complexes (à franges fermées ou ramifiées) conduisant à un nuage de points tridimensionnel S. Ce paragraphe décrit ensuite le traitement de ces sites S par prédiction pour la reconstruction du front d'onde et l'estimation de ses caractéristiques (PTV, RMS etc.) -23- 2778474 La figure 23 représente un réseaux de franges complexe (ici un pointselle). Le traitement décrit ici est le suivant: - une égalisation d'histogramme adaptative, - Une squelettisation des franges, - une recherche des niveaux, - une recherche optimisée des extrema des franges, - un traitement par prédiction. L'égalisation d'histogramme adaptative peut être mise en ouvre de la façon suivante: aux points d'un quadrillage (par exemple 10x10) couvrant la pupille, on calcule un histogramme local sur une fenêtre carrée dont le coté est calculé suivant la distance interfranges locale: par exemple 2 fois cette distance. On trouve un coefficient d'égalisation d'histogramme local. En chaque point de la pupille, on applique alors une correction d'histogramme à partir d'une interpolation des coefficients calculés en chaque point du quadrillage. L'interpolation peut être effectuée par prédiction linéaire avec une fonction d'auto-corrélation choisie arbitrairement, et bien " lissante " c'est à dire de dérivée nulle en 0, puis décroissante. On réalise ensuite une squelettisation des franges d'un type, sombre ou clair, en référence à la figure 24. Il est également possible de traiter les deux types successivement. La recherche des niveaux, c'est à dire des altitudes, de chaque frange est effectuée par un algorithme original dit " par insémination ", comme l'illustre la figure 25. On choisit deux franges voisines. On affecte arbitrairement le niveau 0 à l'une (appelée frange des appuis), et le niveau 1 à l'autre (appelée frange des pivots). Pour chaque point pivot, on prend la tangente à la frange, puis la normale, orientée de l'appui vers le -24- 2778474 pivot. On envoie une " graine de niveau 2 " suivant la normale. Si cette graine atteint un squelette de frange voisine, elle insémine cette frange avec un niveau 2: l'un de ses points possède un " poids " égal à 2, et une " force " inversement proportionnelle à la distance parcourue par la graine. A la fin de l'insémination par le pivot considéré, on fait le bilan de l'insémination de chaque frange. Si l'une d'elles est majoritairement inséminée avec des niveaux 2 (moyenne des " poids " de chaque point multiplié par sa " force "), et minoritairement par des niveaux 0, on décide que son niveau est 2. On propage ensuite l'insémination: chaque frange nouvellement inséminée peut devenir appui ou pivot, et dans ce dernier cas inséminer ses voisines. Cet algorithme est appliqué dans un logiciel de réseaux de franges. On effectue alors une recherche optimisée des extrema des franges, une fois les niveaux reconnus. Cette opération est délicate et cruciale pour la précision du résultat. Comme les franges ont en général un profil asymétrique dans un plan normal vertical, on ajuste à la coupe de l'interférogramme par un tel plan, un polynôme du 3ème degré Avec un traitement par prédiction, on obtient alors un nuage tridimensionnel de points mesurés, que l'on traite par prédiction linéaire pour obtenir la densité de probabilité conditionnelle, et pour l'estimation des caractéristiques du front d'onde, en particulier le PTV, le RMS et le MTF qui est égal à la corrélation normalisée Cor(T) t Cor(0). Une autre application possible est de partir d'un front d'onde déjà reconstruit sur un maillage de nxp valeurs et d'en raffiner l'estimation des -25- 2778474 caractéristiques en prédisant un réseau plus fin de (kn)x(kp) valeurs. Le procédé selon l'invention trouve aussi une application dans l'optimisation de la mesure et de la commande par la méthode de ShackHartmann. En mesure, cette méthode connue d'analyse de surface d'onde [10] utilise une matrice de micro-lentilles. Pour un front d'onde plan, chaque lentille focalise la lumière en son centre Si. Si le front d'onde n'est pas plan, le déplacement du barycentre de la tache focale est réputé représenter la pente locale du front d'onde. En commande, le miroir possède un nombre a d'actionneurs et un nombre c de capteurs (chaque capteur est constitué d'une micro- lentille et d'une sous-pupille du CCD). Les quantités c et a doivent respecter la condition c>a. La relation entre le vecteur U des tensions de commande des actionneurs et le vecteur Y' des pentes mesurées, est linéaire: Y'=L. U (1.15) o L est une matrice rectangulaire dite d'interaction. Les colonnes Lk de L sont les vecteursY de pentes mesurées pour des vecteurs de commande Uk égaux aux vecteurs de la base canonique (1,0,...),(0,1,O,...),(0,0,1,...) etc. La mesure de ces pentes Lk, pour déterminer la matrice L, est effectuée une fois pour toutes. * En utilisation, on mesure un vecteur de pentes Ys et on commande le miroir par un vecteur de tensions U qui est solution des moindres carrés de l'équation (1.15): U=L.Y' =(1.16) o L est l'inverse généralisée de L. Il est à noter qu'en l'absence de prédiction, l'équation (1.15) prend la forme: Us=Ls+.Y) -26- 2778474 Le procédé de reconstruction selon l'invention dans son application à la prédiction des altitudes Z(X) conditionnées par les pentes Z' (S) mesurées, peut être mis en oeuvre pour optimiser l'utilisation de la méthode de Shack-Hartmann en mesure. En particulier pour la reconstruction du front d'onde et pour l'estimation des caractéristiques PTV et RMS du front d'onde et de ses " pentes au pas T ", 8ZT. Pour optimiser l'utilisation de la méthode de Shack- Hartmann en commande des miroirs déformables, on propose de prédire les pentes Y' du front d'onde sur un maillage X plus fin que celui (noté S) des micro-lentilles, en référence à la figure 27. En utilisant les résultats précédents de prise en compte des nonstationnarités, on exprime Y' dans l'équation (1.16) par: Y,= I4 y. (1.17) Dans le terme W, la corrélation Corz doit être estimée en fonction des données Y . Pour une commande rapide, en particulier pour l'optique adaptative, il est possible de n'estimer Corz qu'en temps différé et à de grands intervalles de temps, par exemple toutes les minutes. En effet, chaque trajectoire Z(co) est la réalisation d'une fonction aléatoire Z considérée comme stationnaire, au moins sur des intervalles de temps grands par rapport à la période de commande. Dans L, les colonnes Lk s'expriment alors par L (I4, Lk o Wk ne dépend que de S, X, f,et (CorZ)k. Il est en effet raisonnable d'estimer Corz pour chacune des mesures k, et donc des surfaces d'onde, ce qui ne pose pas de -27- 2778474 problème de temps de calcul puisque L est estimé une fois pour toutes, hors utilisation. Alors en écrivant un partitionnement de L en L+=(L+s, Lx) (1.18) et d'après (1.17) et (1.18): U=L+.Y'=(L+S L+X [)y, on obtient enfin: U = (L+s + Lx. W). Yà (1.19) La commande avec prédiction linéaire ajoute à la matrice de commande classique Ls le terme correctif Lx.W, ce qui ne modifie pas le temps de calcul du vecteur U des tensions de commande du miroir déformable. Le procédé de reconstruction selon l'invention peut avantageusement être appliqué dans les domaines suivants: - la métrologie tridimensionnelle des pièces mécaniques; - la diminution du temps de mesurage en rugosimétrie des pièces mécaniques; l'optimisation de reconstruction de la surface et de l'estimation des paramètres; - l'optimisation de la mesure des surfaces d'onde par interferométrie et l'estimation, de leur PTV, RMS et MTF; l'optimisation de la mesure des surfaces d'onde par méthode de Shack-Hartmann; - l'optimisation de la commande des miroirs déformables dans une chaîne d'asservissement dont le capteur est un dispositif de Shack-Hartmann. Il est à noter que le procédé de reconstruction selon l'invention peut être appliqué avec d'autres technologies de mesure des pentes d'un front d'onde, par exemple l'interférométrie différentielle (ou méthode de Nomarski, en microscopie), ou la méthode des franges de moiré. -28- 2778474 Bien sûr, l'invention n'est pas limitée aux exemples qui viennent d'être décrits et de nombreux aménagements peuvent être apportés à ces exemples sans sortir du cadre de l'invention. -.29 - 2778474 REFERENCES [ 1] CRESSIE N. 1990 The origins of Kriging. Mathematical Geology. Vol.22, No.3. [ 2] DUMAS DE RAULY D. 1968. L'estimation statistique. Gautiers-Villars. Paris. [ 3] FOURGEAUD & FUCHS 1972. Statistique. Dunod. Paris. [4] MATHERON G. 1963. Principles of geostatistics. Economic geology 58, 1246-1266. [ 5] MESTRE M., ABOU-KANDIL H..1993 Linear prediction of signals applied to dimensional metrology of industrial surfaces, Measurement 11, 119-134. Elsevier, Amsterdam. [ 6] MESTRE M., ABOU-KANDIL H., "Measuring the errors of form of industrial surfaces: prediction and optimization", IPrecision Engineering vol. 16 n 4, Butterworth-Heinmann (oct. 1994). [ 7] MESTRE M. 1995 " Méthode des moindres carrés ajustant un modèle déterministe à des données expérimentales avec prédiction linéaire. Estimation, avec prédiction, d'un intervalle de confiance du R.M.S." International Workshop on Mirror Substrate Alternatives, Grasse, Octobre 95. [ 8] ROUBINE E.1979. Introduction à la théorie de la communication. Masson. [ 9] SACKS J. et al. 1989. Design and Analysis of Computer Experiments. Statistical Science Vol.4, n 4, 409-435. [ 10] SOUTHWELL W.H. Wave-front estimation from wave-front slope measurement. J. Opt Soc. Am., Vol 70. No. 8, Aug. 1980. [ 11] YLVISAKER D. 1987. Prediction and design.Ann. Stat. 15, 1-19. -30- 2778474 REVENDICATIONS 1. Procédé pour reconstruire une surface (Y) connue par un ensemble de mesures en des sites (S), mettant en oeuvre une prédiction linéaire en des points (X) situés en dehors des sites (S), caractérisé en ce qu'il comprend une utilisation de mesures de variables géométriques [Y(S), Y'(S)] associées à ladite surface (Y), et en ce qu'il procure une estimation optimale de caractéristiques desdites variables géométriques [Y(S), Y'(S)]. 2. Procédé selon la revendication 1, caractérisé en ce qu'il comprend une utilisation de mesures d'altitudes (Y(S)), et en ce qu'il procure une estimation optimale des caractéristiques géométriques de la surface (S). 3. Procédé selon la revendication 2, mis en ouvre pour la prédiction d'une fonction d'altitude Z(X) conditionnée par des altitudes Z(S) mesurées, caractérisé en ce qu'il comprend en outre un affinage de l'estimation de la fonction d'auto-corrélation et de la fonction de variogramme de la fonction aléatoire Z, cet affinage comprenant une utilisation de la prédiction effectuée à partir d'une première estimation grossière desdites fonctions d'auto-corrélation et de variogramme. 4. Procédé selon la revendication 3, caractérisé en ce qu'il comprend en outre une seconde étape de prédiction. 5. Procédé selon la revendication 1, caractérisé en ce qu'il comprend une utilisation de mesures de pentes (Y' (S)), et en ce qu'il procure une estimation optimale des caractéristiques des pentes (Y'(S)). - 31 - 2778474 6. Procédé selon la revendication 5, mis en oeuvre pour la prédiction des altitudes Z(X) conditionnées par des pentes Z'(S) mesurées, caractérisé en ce qu'on adapte5 aux pentes une solution des moindres carrés généralisée: = (Ifs,.ffltffCYf 7. Procédé selon l'une des revendications 5 ou 6, mis en oeuvre pour la prédiction des pentes Z'(X) conditionnées par des pentes Z'(S) mesurées, caractérisé en ce qu'il comprend en outre une détermination d'un vecteur complet des pentes mesurées et prédites. 8. Procédé selon l'une quelconque des revendications précédentes, appliqué à la rugosimétrie tridimensionnelle de pièces mécaniques, dans lequel on réalise une reconstitution visuelle par tirage aléatoire d'une réalisation de Z(X) avec comme densité la densité de probabilité conditionnelle de (ZxlZs). 9. Procédé selon la revendication 8, caractérisé en ce qu'il comprend les étapes suivantes: - on effectue une simulation avec des points X disposés à l'intérieur d'un premier pavé bordé par deux profils mesurés, puis - on simule dans un second pavé, limitrophe du premier pavé, avec comme points-conditions ceux des profils mesurés qui bordent le second pavé, ainsi que ceux des points X simulés à l'étape précédente et qui bordent aussi le second pavé. -32- 2778474 10. Procédé selon l'une quelconque des revendications 1 à 7, appliqué à la métrologie tridimensionnelle de pièces mécaniques. 11. Procédé selon l'une quelconque des revendications précédentes, mis en oeuvre dans un système d'analyse de surface d'onde par interférométrie, ce système comprenant une analyse de réseaux de franges d'interférences conduisant à un nuage de points tridimensionnel S, un traitement de ces points S par prédiction pour la reconstruction d'un front d'onde et une estimation des caractéristiques de ce front d'onde, caractérisé en ce que le traitement comprend: - une égalisation d'histogramme adaptative, - Une squelettisation des franges, - une recherche des niveaux, - une recherche optimisée des extrema des franges, et - un traitement par prédiction. 12. Procédé selon la revendication 11, caractérisé en ce que l'égalisation d'histogramme adaptative comprend les étapes suivantes: - on calcule un histogramme. local sur une fenêtre carrée dont le coté est calculé suivant la distance interfranges locale, aux points d'un quadrillage couvrant une pupille, - on en déduit un coefficient d'égalisation d'histogramme local, - on applique alors une correction d'histogramme à partir d'une interpolation des coefficients calculés en chaque point du quadrillage. 13. Procédé selon la revendication 12, caractérisé en ce que l'interpolation est effectuée par prédiction linéaire avec une fonction d'auto-corrélation choisie arbitrairement. 33 2778474 14. Procédé selon l'une quelconque des revendications 11 à 13, caractérisé en ce qu'on réalise une squelettisation des franges d'un type, sombre ou clair. 15. Procédé selon l'une quelconque des revendications 11 à 14, caractérisé en ce que la recherche des niveaux de chaque frange est effectuée par un algorithme par insémination. 16. Procédé selon la revendication 15, caractérisé en ce que l'algorithme par insémination comprend les étapes suivantes: - on choisit deux franges voisines, - on affecte arbitrairement un niveau 0 à l'une desdites franges appelée frange des appuis, et un niveau 1 à l'autre frange appelée frange des pivots, pour chaque point pivot, on prend la tangente à la frange, puis la normale, orientée de l'appui vers le pivot, - on envoie une graine de niveau 2 suivant la normale, - si cette graine atteint 'un squelette de frange voisine, elle insémine cette frange avec un niveau 2. 17. Procédé selon la revendication 16, caractérisé en ce qu'à la fin de l'insémination par le pivot considéré, on réalise un bilan de l'insémination de chaque frange, de sorte que, si l'une desdites franges est majoritairement inséminée avec des niveaux 2, et minoritairement par des niveaux O, on lui attribue un niveau 2. 18. Procédé selon la revendication 17, caractérisé en ce qu'il comprend en outre une opération de propagation -34 -2778474 de l'insémination, dans laquelle chaque frange nouvellement inséminée peut devenir appui ou pivot. 19. Procédé selon la revendication 18, caractérisé en ce qu'il comprend en outre, après reconnaissance des niveaux, - une recherche optimisée des extrema des franges, - un traitement par prédiction à l'issue duquel on obtient un nuage tridimensionnel de points mesurés, - un traitement par prédiction linéaire des points mesurés pour obtenir une densité de probabilité conditionnelle et pour estimer des caractéristiques du front d'onde. 20. Procédé selon la revendication 11, caractérisé en ce qu'on part d'un front d'onde déjà reconstruit sur un maillage de nxp valeurs et en ce qu'on raffine l'estimation des caractéristiques en prédisant un réseau plus fin de (kn)x(kp) valeurs. 21. Procédé selon l'une quelconque des revendications 1 à 7, appliqué à la prédiction des altitudes Z(X) conditionnées par les pentes Z' (S) mesurées, caractérisé en ce que ce procédé est appliqué à une optimisation d'une opération de mesure réalisée par la méthode de Shack-Hartmann. 22. Procédé selon la revendication 21, caractérisé en ce qu'il est en outre appliqué à une optimisation d'une opération de commande réalisée par la méthode de Shack- Hartmann. 23. Procédé selon la revendication 22, mis en oeuvre pour optimiser l'utilisation de la méthode de Shack- Hartmann en commande des miroirs déformables, cette -35- 2778474 méthode utilisant une matrice de micro-lentilles, caractérisé en qu'il comprend une prédiction des pentes (Y') du front d'onde (Y) sur un maillage plus fin que celui des micro-lentilles.