FR 2815133 A1 20020412 FR 0012995 A 20001011 PROCEDE DE TRAITEMENT EN VUE DE LA LOCALISATION D'UNE SOURCE GENERATRICE D'ONDES, NOTAMMENT D'ONDES SISMIQUES DOMAINE GENERAL DE l'INVENTION La présente invention concerne un procédé de traitement en vue de la localisation, dans le milieu dans lequel elle se propage, d'une source génératrice d'ondes (tel qu'un événement sismique dans son milieu géologique), observée par au moins un récepteur multi-axial (tel qu'un réseau de capteurs sismiques permanent ou temporaire, déployé en surface ou en profondeur, etc.). En sismologie, la localisation des tremblements de terre à partir des sismogrammes enregistrés sur un réseau sismique est omniprésente. Cependant, dans la majorité des cas, cette localisation est estimée en supposant que le milieu de propagation entre la source et les récepteurs sismiques est relativement simple : milieu à vitesse constante, milieu latéralement homogène mais à couches horizontales.... Dans de tels milieux, la description analytique de la propagation des ondes est relativement simple et permet de localiser le tremblement de terre à l'aide de méthodes inverses qui sont appliquées directement aux équations de propagation ou bien qui utilisent des processus d'itération. En milieu de propagation complexe, ou 3D, il n'est actuellement pas possible d'inverser les équations de propagation des ondes sismiques. Il existe des méthodes itératives, mais d'une part celles-ci ne permettent pas d'évaluer la distribution des localisations possibles, puisqu'elles ne donnent que la meilleure localisation, et d'autre part, elles peuvent être instables. Un but particulier de l'invention est de proposer un procédé de traitement qui permette d'évaluer avec une grande fiabilité une distribution de probabilités de localisation de l'épicentre d'un séisme, ce traitement pouvant toutefois être étendu de façon plus générale à la distribution de probabilités de toute source génératrice d'ondes. ETAT DE LA TECHNIQUE Un exemple de signal enregistré sur une station à trois composantes orthogonales, suite à la propagation des ondes sismiques engendrées par un micro-séisme, est illustré sur la figure 1. Dans tout le présent texte, les terminologies de récepteur ou de station sont utilisées indifféremment pour désigner des capteurs sismiques à trois axes. Par ailleurs, un événement sismique ou microsismique ou bien un tremblement de terre ou encore un séisme sont ici synonymes. Il a déjà été proposé des traitements permettant d'utiliser l'information des sismogrammes enregistrés par les différents récepteurs d'un réseau afin d'obtenir une distribution de la probabilité de répartition de l'épicentre sismique. On pourra notamment à cet égard se référer aux articles et publications suivants : Lomax, A., 2000. http : //iapetus. unice. fr : 80/-lomax/nlloc. Lomax, A., Virieux, J., Volant, P. & Berge, C., 2000. Probabilistic earthquake location in 3D and layered models. Introduction of a Metropolis-Gibbs method and comparison with linear locations, in Advances in Seismic Event Location. Kluwer, Amsterdam. Accepted, in press. Nelson, G. D. & Vidale, J. E., 1990. Earthquake locations by 3D finite-difference travel times, Bull. Seis. Soc. Am., 80 (2), 395-410. Tarantola, A. & Valette, B., 1982. Inverse problems = quest for information, J. Geophys., 50,159-170. Tarantola, A., 1987. Inverse problem theory. Elsevier, New-York. Wittlinger, G., Herquel, G. & Nakache T., 1993. Earthquake location in strongly heterogeneous media, Geophys. J. Int., 115,759-777. Les traitements présentés dans ces différents articles et publications mettent en oeuvre une détermination à partir des temps d'arrivée des ondes aux récepteurs. Ces temps d'arrivée pointés sont généralement ceux qui correspondent à l'onde première compressive P et à l'onde seconde de cisaillement S observées sur les sismogrammes enregistrés par les récepteurs sismiques, phases les plus facilement identifiables. Sur la figure 1, ces pointés sont matérialisés par les deux traits gras. Il est en effet possible, suivant l'approche de Tarantola et Valette (1982) qui utilisent une norme L2 et se placent sous le critère des moindres carrés, d'écrire que la fonction densité de probabilité (fdp) calculée a posteriori sur la localisation spatiale X (correspondant à trois coordonnées x, y, z) de l'événement, estimée à partir de l'ensemble des temps d'arrivée des phases sismiques est : a, (X) = K p, (X) exp- (t.-h' (x))' (C. +C, " (t..-h' (X)) l) Equation 1, 2 11 avec : a K : facteur de normalisation introduit lors de l'intégration sur le temps. (Il faut en effet que l'on +00 puisse avoir J'ot (X, t) dt = 1). oxo . pt (X) : fdp liée à l'information dont on dispose a priori sur les données et la modélisation de la propagation, cette densité pouvant dépendre de la position X de la source. 'to'est un vecteur qui regroupe les temps d'arrivée de toutes les phases identifiées à toutes les stations et qui est défini par les composantes La composante top étant égale pour une onde sismique donnée et une station i donnée au temps observé à la station i moins la moyenne pondérée des temps observés à toutes les stations, les coefficients pj de pondération étant définis par p, = -, (où PJi désigne la composante 1 de la jème ligne et de la ième colonne de la matrice P) 'h' (X) étant quant à lui un vecteur dont les composantes h sont définies par h i (X) étant égal pour une onde sismique donnée et une station donnée, au temps de propagation jusqu'au point X modélisé depuis la station i moins la moyenne pondérée des temps modélisés à toutes les stations. Ces temps d'arrivée théoriques peuvent notamment être déterminés, à partir du modèle de vitesse de propagation de l'onde considérée, par calcul des paramètres de rais tirés à partir de sources réparties sur les noeuds d'une grille 3D qui échantillonne une zone cible où l'on estime que l'épicentre sismique est susceptible de se trouver ; à cet égard, la figure 2 illustre un rai de propagation P entre une source S et un récepteur R, ainsi qu'une grille G d'échantillonnage d'une zone cible définie par une pluralité de noeuds N. Cette méthode d'estimation des temps théoriques d'arrivée des ondes au (x) récepteur (s) est ce qu'on appelle une méthode par tracé de rais . . (Ct +CTr'=pl où : - Ct : matrice de covariance calculée sur les temps d'arrivée des phases sismiques observés au niveau des différentes stations du réseau ; - CT : matrice de covariance calculée sur les temps d'arrivée théoriques déterminés en fonction d'un modèle de vitesse de propagation dont on dispose. Il est ainsi possible, en donnant à X des valeurs discrètes correspondant aux coordonnées de noeuds d'une grille par laquelle on échantillonne une zone cible et moyennant la connaissance de temps de trajets estimés ou modélisés entre lesdits noeuds et les récepteurs du réseau, de calculer une répartition de probabilités de localisation ot (X) pour les noeuds X de ladite grille. Le temps T correspondant à l'origine de l'événement sismique peut également être estimée a posteriori par la formule : avec, pour une phase sismique donnée : N : nombre de pointés disponibles sur le réseau, tel : pointé effectué à la station i, hi (Xi) : temps de propagation entre la station i et la source placée en Xi, modélisé (par tracé de rais, par exemple). PRESENTATION DE L'INVENTION L'invention propose quant à elle une autre technique d'évaluation de distributions de probabilités de la localisation d'une source d'onde (s). Cette technique permet en elle-même une localisation particulièrement précise de cette source. Elle est en outre avantageusement utilisée en complément de la technique de détermination de la distribution de probabilités de localisation obtenue à partir des temps d'arrivée des ondes au (x) récepteur (s). Ainsi l'invention propose un procédé probabiliste de localisation de la source d'ondes qui prend en compte la polarisation des ondes observées au niveau des récepteurs pour en déduire une densité de probabilité de localisation de cette source dans une zone cible. Par polarisation d'une onde on entend la mesure, en un point de l'espace, de la répartition de l'énergie (ou de l'amplitude) de cette onde dans l'espace. Par vecteur polarisation, on entend le vecteur pointant dans la direction (et non le sens) du maximum d'énergie. On notera à cet égard que la polarisation des ondes apporte une indication sur la direction dans laquelle se trouve la source, vue du récepteur, tandis que la mesure des temps d'arrivée des ondes fournit une indication sur l'éloignement de la source par rapport au récepteur. Plus précisément, l'invention propose un procédé de détermination d'une distribution de probabilités de localisation d'une source d'ondes captées par au moins un récepteur multi-axial, selon lequel : - on détermine à partir des signaux relevés par ledit récepteur la polarisation des ondes au niveau de celuici, - on détermine pour une pluralité de points d'une zone cible où la source est susceptible de se trouver une polarisation théorique attendue audit récepteur pour le cas où la source se serait trouvée au point de la zone cible considérée, - on calcule pour les différents points considérés dans la zone cible une probabilité de localisation qui est fonction de la polarisation attendue au récepteur et de la polarisation mesurée. En conséquence, il est indispensable que la polarisation observée puisse être modélisée, cette modélisation étant elle-même dépendante de la position de la source. Notamment, selon un mode de mise en oeuvre préféré, on calcule pour les différents points considérés dans la zone cible une probabilité de localisation qui est fonction, d'une part, des composantes d'une matrice de covariance ou d'auto-corrélation calculée pour une onde donnée à partir des signaux enregistrés au niveau du récepteur et, d'autre part, des composantes du vecteur polarisation théorique attendu audit récepteur, pour cette même onde, et pour le cas où la source se serait trouvée au point de la zone cible considérée. En particulier, dans les cas où le signal généré par la source est observé en une pluralité de récepteurs et/ou plusieurs ondes et donc polarisations sont observées (et modélisées) pour un même récepteur, on calcule pour les différents points de la zone cible considérés une probabilité-conjointe de localisation qui est le produit des différentes probabilités de localisation pour les différents récepteurs et/ou les différentes polarisations qui y sont observées. Le procédé de détermination de la distribution des probabilités à partir des données de polarisation est notamment avantageusement utilisé conjointement à une détermination de la distribution des probabilités à partir des données de temps d'arrivée, ce qui permet d'exploiter les redondances contenues dans les informations sur les temps d'arrivée et sur les polarisations. En multipliant les probabilités obtenues à partir des unes et des autres de ces données, on obtient une distribution de probabilité-conjointe qui permet d'appréhender de façon plus fine la localisation de la source. Ainsi, selon un mode de mise en oeuvre avantageux de l'invention : - on détermine à partir des signaux relevés par un ou plusieurs récepteurs le ou les temps d'arrivée des ondes au niveau de celui-ci, - on détermine pour les différents points constituant la zone cible un ou des temps d'arrivée théoriques attendus audit récepteur pour le cas où la source se serait trouvée au point considéré de la zone cible, - on détermine pour tous les points de la zone cible une probabilité de localisation qui est le produit, d'une part, d'une probabilité de localisation fonction des polarisations attendues et mesurées au (x) récepteur (s) et, d'autre part, d'une probabilité de localisation fonction des temps d'arrivée attendus et mesuré au (x) récepteur (s). On notera que le procédé proposé par l'invention est particulièrement avantageux dans le cas où l'on dispose d'un réseau comprenant peu de récepteurs, la polarisation permettant, malgré le faible nombre de récepteurs, de contraindre fortement la position de la source et donc d'obtenir une certaine précision sur cette position. Egalement, on ne possède pas toujours un modèle de propagation de chaque onde dans la zone cible. Par exemple pour l'application sismologique, si on ne dispose que d'un modèle de propagation pour les ondes P, on ne peut pas exploiter les temps d'arrivée des ondes S observées au niveau des récepteurs. Dans ce cas, le procédé de l'invention est particulièrement avantageux puisqu'il permet de préciser les résultats grâce aux données de polarisation des ondes P. PRESENTATION DES FIGURES D'autres caractéristiques et avantages de l'invention ressortiront encore de la description qui suit, laquelle est purement illustrative et non limitative et doit être lue en regard des dessins annexés sur lesquels : la figure 1, déjà discutée, est un exemple de sismogrammes enregistrés sur une station à trois composantes orthogonales : Z, X et Y ; la figure 2, également déjà discutée, est une représentation schématique illustrant un rai de propagation (P) déterminé par modélisation entre un point source (S) et un point récepteur (R) et une grille d'échantillonnage d'une zone cible définie par une pluralité de noeuds (N) ; - la figure 3 est une illustration géométrique de la matrice de covariance d'un signal trois composantes dans son repère propre ; la figure 4 est une représentation schématique illustrant une détermination d'incertitude de localisation ; - la figure 5 est une représentation schématique de la localisation d'une source sous forme de lignes d'isoprobabilités. DESCRIPTION D'UN OU PLUSIEURS MODES DE MISE EN OEUVRE DE L'INVENTION 1. Généralités 1.1 RAPPELS : DISTRIBUTION GAUSSIENNE On rappelle qu'une distribution gaussienne associée à une variable X quelconque s'écrit : g (X) = A exp (x-X)'C-' (x-x) 1j 2 Il Une telle distribution peut s'interpréter comme une fonction densité de probabilité (fdp) si son intégrale sur l'espace supporté par la variable X est égale à 1. La constante A est une constante de normalisation. Le terme de droite dans la formule ci-dessus, la constante A mise à part, est appelé fonction caractéristique. L'argument de l'exponentielle est appelé fonction coût ou fonction écart. 1.2. POLARISATION D'ONDE L'invention met en oeuvre un traitement à partir de la polarisation d'ondes. L'invention trouve en particulier avantageusement application dans le cas d'ondes élastiques et notamment d'ondes sismiques ; elle pourrait toutefois également touver à s'appliquer dans le cas d'ondes électromagnétiques, ou d'ondes qui peuvent être décrites par la théorie de l'optique. Dans l'application sismique, la polarisation de la phase P (onde première compressive) est la plus accessible sur des données réelles. Par ailleurs, l'onde P, toujours tangente au rai sismique qui relie la source au point où elle est observée, est, de ce fait, liée à la direction dans laquelle se trouve la source par rapport au point où cette onde est observée. La polarisation de la phase S (onde seconde de cisaillement) est plus souvent sujette (car plus sensible) à des effets de propagation mal modélisés ou mal estimés comme l'anisotropie, ce qui la rend plus difficilement exploitable. Si ce n'est pas le cas, elle est polarisée dans un plan perpendiculaire au rai sismique. 1.2. 1 Matrice de covariance Pour déterminer la polarisation d'une onde, ici sismique, on utilise avantageusement la méthode décrite dans Kanasewich, E. R., (1973). Time sequence analysis in geophysics. The University of Alberta Press, Edmonton. basée sur l'analyse de la matrice de covariance. Cette matrice est calculée à partir du signal enregistré sur les trois composantes orthogonales du sismomètre et pour une fenêtre temporelle choisie, cette fenêtre étant de préférence une fenêtre centrée sur les ondes qui nous intéressent et en particulier l'onde P (la formulation étant néanmoins quelle que soit la position de la fenêtre). Cette méthode pré-suppose une distribution gaussienne de chacun des échantillons du signal superposé à du bruit blanc, le calcul étant réalisé en se plaçant (comme dans le cas des données de type temps d'arrivée) dans le domaine de la norme L2. On se place dans la suite dans le cas d'un repère cartésien orthonormé (x = Est, y = Nord, z = Verticalhaut). Soi t Xij l'échantillon j du signal sur la composante i. La moyenne (ou espérance mathématique) de cette composante sur la fenêtre temporelle composée de n échantillons est : La matrice de covariance (ou d'auto-corrélation) qui mesure la dispersion du signal à moyenne nulle dans l'espace est l'auto-corrélation du signal centré, c'est- à-dire : r- ou encore : 1 n r--J-Vy yVy y) n 1 kil Cette matrice de covariance est symétrique, réelle et définie positive. Il existe donc une base orthonormée directe dans laquelle C est diagonale. Cette base, appelée base propre, est définie à l'aide des vecteurs propres de C. Ainsi, CVi = Vi, Vi vecteur propre de C, de valeur propre Ài. Physiquement, dans la base propre, chacune des valeurs propres est proportionnelle à l'énergie du signal dans la direction du vecteur propre qui lui est associé. Pour une onde première compressive (P), cette énergie est maximale dans la direction d'arrivée du rai à la station et devrait correspondre à la direction du vecteur propre associé à la valeur propre maximale. Du point de vue géométrique, la matrice de covariance prend, dans la base propre, l'allure d'un ellipsoïde. La figure 2 est une représentation de la matrice de covariance dans son repère propre. Les vecteurs propres Vi, V2 et V3 de la matrice de covariance définissent une nouvelle base orthonormée (dans le repère géographique). L'énergie du signal étudié est distribuée sur ces trois axes et est proportionnelle aux valeurs propres correspondantes. L'analyse des valeurs propres les unes par rapport aux autres permet d'étudier la forme de l'ellipsoïde et par analogie la distribution de l'énergie dans l'espace. En supposant que À : JI, on introduit trois nouveaux paramètres : chacun de ces paramètres l, p et s étant compris dans l'intervalle [0, 1], ces différents paramètres vérifiant 1+p+s=1. Il est important de noter que : si 2 = 3 == 0, alors-Z = 1 et p = s = 0, le signal est polarisé linéairement dans la direction Vi, Si X 2 et À3 = 0, alors 1 = s = 0 et p = 1, le signal est polarisé dans le plan orthogonal à V3, Si k 2 = À3, alors 1 = p = 0 et s = 1, le signal est isotrope et n'est donc pas du tout polarisé. On reviendra sur ces propriétés dans la suite du texte. On notera que le fait de savoir a priori que le signal est polarisé dans un plan permet l'utilisation de récepteurs bi-axiaux (et non tri-axiaux) avec les deux composantes orthogonales placées dans le plan de polarisation a priori. La formulation ci-dessus est alors restreinte à un espace à deux dimensions uniquement. Pour rester dans le cas le plus général, nous continuerons à illustrer nos propos avec des récepteurs tri-axiaux. Il convient toutefois de souligner que la notion de polarisation du signal sur une fenêtre temporelle est totalement incluse dans la matrice de covariance du signal calculée sur cette même fenêtre. 2. Exemple de calcul d'une densité de probabilité de localisation 2.1 DENSITE DE PROBABILITE POLARISATION On se place en une station trois composantes donnée. Pour un couple station/noeud de la grille cible (placé en X), on définit la densité de probabilité a posteriori sur la localisation spatiale de l'événement estimée à partir des données polarisation sur la station comme : (X) = (X) exp--p, (x))' (C. +Cj' (P.-P, (x)) 1j Equation 2, 2 avec : Po : le vecteur polarisation observé au récepteur considéré, Pm (X) : le vecteur polarisation modélisé qui devrait être observé au récepteur et qui est par exemple déterminé par tracé de rais, à partir d'un modèle de vitesse dont on dispose, ce vecteur dépendant de la position X de la source, Co : la matrice de covariance sur la polarisation observée, Cm : la matrice de covariance sur la polarisation modélisée. A : coefficient de normalisation qui permet de différencier la densité de probabilité de la probabilité elle-même. Ce coefficient est tel que la somme des probabilités sur l'espace supporté par le vecteur polarisation soit égale à 1. Comme le module de Pm vaut 1, il faut intégrer la fonction caractéristique de la probabilité sur la sphère de rayon 1. A doit donc être tel que : f u, (X, P) dP = 1 Sphere pp (X) : fdp liée à l'information a priori sur les données et le modèle, cette densité de probabilité étant susceptible de dépendre de la position X de la source ; par exemple elle peut-être égale à 1 dans la zone cible et à 0 sinon. Cette équation 2 appelle plusieurs remarques. a) On note que cette équation, le terme en pp (X) mis de côté, représente bien une distribution gaussienne, cf. 1.1. b) Par hypothèse (cf. ≈1.2. 1), la polarisation observée est supposée centrée, donc à moyenne nulle, et sa dispersion entièrement décrite par la matrice de covariance du signal. Par suite, dans l'équation 2, le vecteur Po = O. c) Par ailleurs, on peut supposer que le vecteur polarisation modélisé est exact et éliminer dans l'équation la matrice Cm. Par conséquent, l'Equation 2 se simplifie et devient : U, (X) = A p, (X) exp- PJ (X) C.' P, (X))] ! 1-2 m Il Pour éliminer tout effet d'amplitude dans les données observées, c'est-à-dire dans la matrice de covariance associée à l'onde considérée, et pouvoir ainsi traiter de la même façon toute donnée polarisation, on normalise la matrice de covariance par la plus grande valeur propre. L'équation précédente devient alors : ap (X) = (x) expf- (Pj (X) CP, (X)) j ! Equation 3, 2 1 avec C=-cJL la matrice de covariance normalisée de la Âl il, polarisation observée au récepteur. Dans la base définie par les vecteurs propres de la matrice de covariance des données observées, l'Equation 3 s'écrit encore : fi (X) exp P22 (X) + 2 (X)) ] Equation 4, 2 -12 -13 avec : Pj, P2 et P3 : les composantes du vecteur Pm exprimées dans la base propre de la polarisation observée, la composante Pi est la projection du vecteur Pm sur le vecteur Vi définissant la base propre (i = 1, 2 ou 3) ; À1, À2 et k3 : les valeurs propres de la matrice de covariance calculée sur les données observées. Ainsi que précédemment indiqué, les valeurs propres décrivent la forme de l'ellipsoïde de polarisation et donc la qualité de la polarisation. Il est donc physiquement justifié que les valeurs propres apparaissent dans la densité de probabilité, ce que souligne l'Equation 4. D'après cette équation, on voit et on peut vérifier que, pu (X) mis à part : si i = 2 = À3, le signal n'est pas polarisé (ou est isotrope) et comme M=lj, alors ap (X) est constante, si k 2 et À3 = 0, le signal est polarisé selon un plan perpendiculaire à V3 et ap (X) = 0 sauf sur ce plan où elle est constante (P3 = 0), si k2 = 13 = 0, le signal est polarisé selon Vi et ap (X) = 0 sauf le long de cette direction (P2 = P3 = 0) où elle est constante. Il est ainsi possible, connaissant une polarisation mesurée sur un récepteur et connaissant la polarisation attendue sur ce récepteur selon la position de la source dans une zone cible donnée, de calculer une probabilité de localisation a posteriori de l'événement sismique dans cette zone cible. Bien entendu, pour un séisme donné, plusieurs vecteurs polarisations peuvent être obtenus par sismogramme mais aussi par récepteur. En conséquence, si l'on suppose que ces vecteurs sont indépendants les uns des autres, alors la fdp globale a posteriori associée aux polarisations est : Equation 5 avec M le nombre de polarisations disponibles sur le réseau. 2.2 CALCUL D'UNE DENSITE DE PROBABILITE CONJOINTE La fdp-globale de localisation ainsi calculée à partir des polarisations mesurées peut être combinée avec la fpd-globale de localisation calculée à partir des temps d'arrivée, dans une fdp-conjointe. Ce qu'on appelle une fdp-individuelle est une fonction densité de probabilité a posteriori associée à une donnée particulière (une polarisation, un temps d'arrivée). Une fdp-globale est déduite de plusieurs fdp-individuelles pour un type de donnée (polarisation, temps d'arrivée). Enfin, une fdp-conjointe est déduite de plusieurs fdpglobales. Si les données temps ne sont pas utilisées, la fdp-conjointe est équivalente à la fdp-globale des polarisations. Cette fdp-conjointe est tout simplement le produit de toutes les fdp-individuelles construites à partir des données observées (et modélisées) de type temps d'arrivée et polarisations disponibles, cette simple multiplication étant justifiée par le fait que toutes les données sont considérées comme indépendantes les unes des autres. La fdp-conjointe a posteriori s'écrit : Equation 6, avec : B : constante de normalisation telle que la somme sur l'espace supporté par les variables de la fdp-conjointe soit égale à 1. S'il n'y a pas d'erreur sur les temps de propagation modélisés par le tracé de rais, CT disparaît de l'équation. M : nombre de polarisations disponibles sur le réseau. Rappelons que les vecteurs to ou h'sont de dimension N, le nombre de pointés disponibles sur le réseau sismique. Donc, pour le séisme, N+M données sont déduites des sismogrammes disponibles sur l'ensemble du réseau. La meilleure localisation est le point de l'espace pour lequel cr (X) est maximale. Remarques a) Il est important de noter que le produit de distributions gaussiennes ne donne pas nécessairement une distribution gaussienne. b) Dans ce qui précède, on s'est limité aux fdpindividuelles calculées à partir des données de type temps et celles de type polarisation. Si une fdp a posteriori avait été définie pour les données de type amplitude des ondes, il aurait été possible d'en tenir compte également. Cependant, dans ce dernier cas, les amplitudes des phases sismiques ne sont pas indépendantes des temps de propagation, et donc des temps d'arrivée des phases. Par conséquent, il serait nécessaire d'introduire des corrélations entre ces deux types de données et la fdp-conjointe ne serait plus le simple produit des fdpindividuelles. 3. Estimation des incertitudes Plusieurs traitements sont possibles pour estimer les incertitudes sur le résultat de la localisation et pour en tenir compte dans la distribution des probabilités. 3.1 INCERTITUDES DEDUITES DES ERREURS SUR LES DONNEES On se place ici du point de vue des données acquises. Pour chacune d'elles, on décide de calculer les fonctions caractéristiques des fdp-individuelles uniquement aux noeuds de la grille cible pour lesquels la variable de la fdp temps ou polarisation est comprise entre plus ou moins un écart type autour de sa valeur espérée. Cela revient à ne garder que les noeuds pour lesquels la fonction coût de la fdp-individuelle est inférieure ou égale à 1. Par suite, uniquement certains noeuds de la grille vont satisfaire cette condition, les autres seront supposés avoir une fdp nulle. a) Ainsi, pour les temps d'arrivée observés, on détermine les points X de la zone cible qui vérifient : (t.-h' (x)) -h' (x)) . Dans le cas où les temps observés et modélisés sont tous indépendants les uns des autres-ce qui n'est pas nécessairement le cas (e. g. lorsque l'on traite l'onde P et l'onde S au niveau d'un même récepteur)-alors la matrice P est diagonale et il suffit pour chacun des temps de déterminer les points X de la zone cible qui vérifient (cf. Equation 1) : b) Pour les polarisations, on peut supposer, comme pour obtenir u,, p (X), que celles-ci sont indépendantes les unes des autres, c'est-à-dire d'un récepteur à l'autre. En conséquence, il faut déterminer les points X de la zone cible qui vérifient (cf. Equation 4) : Â Â P -2 et -en même temps, Âl Â, rappelons que Pl, P2 et P3 sont les composantes du vecteur Pm exprimées dans la base propre de la polarisation observée. c) Par suite, la fdp-conjointe ce (X) (cf. Equation 6) ne sera définie que sur les points qui vérifient l'ensemble des contraintes ci-dessus, pour l'ensemble des récepteurs, des temps d'arrivée et des polarisations mesurées sur ceux-ci ; autrement dit, les points qui résultent de l'intersection de tous les sous ensembles définis par chaque fdp-individuelle. C'est ce qu'illustre la figure 4 sur laquelle on a représenté une zone d'incertitude de localisation autour d'une meilleure localisation ML, cette zone l étant délimitée par les droites EP d'écart type sur les données de polarisation pour deux récepteurs RI et R2, ainsi que par les courbes ET d'écart type sur les données temps pour ces mêmes récepteurs R1 et R2. La distribution de ces points autour de la meilleure solution (celle qui présente la fdp-conjointe maximum) constitue l'erreur de localisation déduite de l'incertitude sur les données. L'analyse de cette distribution permet de quantifier les incertitudes sur la localisation. En supposant que cette distribution peut-être approximée par une distribution gaussienne dans l'espace physique, alors le calcul de la matrice de covariance sur les points qui la constitue permet d'obtenir les vecteurs propres et valeurs propres qui décrivent la distribution. Exactement comme dans le cas du calcul de polarisation (cf. ≈1.2. 1). Pour estimer dz écart type associé à la meilleure localisation et qui englobe tous les points de la distribution, nous projetons tous ces points sur les vecteurs propres et gardons, pour chacun des axes, la distance maximale qui sépare deux points projetés. Cette distance équivaut à deux écarts types dans la direction des vecteurs propres. (Notons que l'amplitude des valeurs propres ne correspond pas à l'écart type de la distribution dans la direction de son vecteur propre associé puisque les points autour de la meilleure localisation appartiennent déjà au domaine de 1 écart type autour de la valeur moyenne.) Cette première définition de l'incertitude sur la localisation permet de s'affranchir du calcul des constantes de normalisation des fdp-individuelles ou des fdp-globales et donc de travailler directement avec leurs fonctions caractéristiques ou bien leurs fonctions coût. Remarques - Comme nous l'avons précédemment fait remarquer, le produit de distributions gaussiennes ne donne pas nécessairement une distribution gaussienne. Par suite, l'étude de la distribution des localisations par la méthode de covariance (cf. paragraphe précédent) n'a de sens physique que si cette distribution est proche d'une distribution gaussienne. En toute rigueur, il serait nécessaire de vérifier que la distribution étudiée est gaussienne ou presque, à l'aide de tests statistiques de Y2 type ou Kolmogorov. Une condition nécessaire mais insuffisante est que le centre de gravité de la distribution (moyenne spatiale) corresponde à la meilleure localisation. - Il peut arriver que l'intersection de tous les groupes de points entre eux soit vide. Deux raisons à cela peuvent être proposées : 'la zone cible n'est pas assez grande et la meilleure localisation se trouve en dehors de celle-ci, 'les données observées sont tellement incohérentes entre elles que la prise en compte d'un écart type n'est pas suffisant pour obtenir une intersection. Dans ce dernier cas, le remède consiste à élargir le critère de sélection sur l'écart entre les données observées et les données modélisées et donc à considérer des divergences supérieures à un écart type. 3.2 INCERTITUDES A POSTERIORI Les incertitudes a posteriori sur la localisation sont déduites de la distribution finale de la fdp-conjointe, or (X), sans sélection a priori sur la valeur des fonctions coûts aux noeuds de la grille. Dans ce cas, contrairement à l'estimation précédente, il est nécessaire de correctement estimer les fdp-globales temps et polarisation pour pouvoir les combiner (si nécessaire) et ensuite pouvoir travailler sur la distribution de la valeur absolue (par opposition à relative) de la fdpconjointe sur la zone cible. On sait comment calculer la constante de normalisation K de la fdp-temps ut (X), après avoir supposé que #t (x) est constante (Tarantola, 1987), si bien que l'Equation 1 devient : 1 0", (x) = ( (2) det (C, + C,)) xp [--h' (X)) (C.-. C,)-' (-h' (X))] ! 2 il avec N la dimension du vecteur te'ou h'. Cependant, pour la fdp-polarisation un problème se pose. Ce qui importe est de calculer correctement la fdp, toutes polarisations confondues, autrement dit o-p (X) (cf. Equation 5). Cela passe donc par le calcul de la constante A de chaque fdp-polarisation ap (X). Or, l'intégrale sur la sphère de rayon 1 ne se calcule pas de façon triviale. Afin de contourner cette difficulté et de traiter de la même façon (si nécessaire) les fdp temps et polarisation, on introduit une fdp a priori sur les données ou le modèle, c'est-à-dire que/ (X) et pu (X) ne sont plus constantes et on passe d'un milieu (ou espace physique) continu à un milieu discrétisé, celui de la grille de la zone cible par exemple. Plus précisement, pour estimer la constante de normalisation de chaque fdp (temps ou/et polarisation) numériquement, on somme sur un ensemble (à définir) de noeuds de la grille cible la valeur de la fonction caractéristique. Ce passage d'une intégrale à une somme discrète introduit un effet d'échantillonnage de l'espace physique réel par la grille cible. La conséquence en est une surestimation de la fdp par rapport à la valeur qu'elle aurait dans le domaine continu. L'ensemble des noeuds sur lequel sommer est définit par la fdp a priori sur les données ou le modèle, fdp a priori qui est choisit de façon arbitraire. Plusieurs possibilités peuvent être proposées : 1) On peut notamment supposer que en dehors de la zone cible, il n'y a plus de localisation possible, que l'on utilise les temps d'arrivée au (x) récepteur (s), ou les polarisations. Cette proposition de fdp a priori signifie encore qu'on a toutes les chances d'avoir la localisation dans la zone cible. Ainsi, On note que si chaque pp (X) vérifie cela, alors la fdp a priori associée à o-çp (X), vérifie aussi cela. Pour normaliser cgp (X) et at (X), il faut donc intégrer leurs fonctions caractéristiques sur l'ensemble de la grille cible. Ceci est fait numériquement en sommant la valeur de la fonction caractéristique sur l'ensemble des noeuds qui constituent la grille, puis en divisant chacune des valeurs aux noeuds par cette somme. 2) D'autres propositions peuvent être faites. On peut par exemple supposer qu'il n'y a plus aucune chance de trouver la localisation après tant de mètres de propagation de l'onde. Dans ce cas, ce qui borne la fdp a priori pp (X), et a fortiori Pgp (X) ou/et #t (X) est une certaine distance curviligne constante, distance mesurée à partir du (des) récepteur (s) le long du trajet de l'onde considérée. Ainsi, L'indice c indique qu'il s'agit d'une distance curviligne. Une autre proposition serait de dire qu'il n'y a plus aucune chance de trouver la localisation au-delà de tant de temps de propagation de l'onde. Dans ce cas, ce qui borne la fdp a priori pp (X) ou/et pt (X) est un certain temps de propagation constant. Pour appliquer ces deux propositions, il est nécessaire de fournir la valeur de la distance parcourue le long du rai par l'onde, ou bien le temps de propagation pour l'onde considérée (nécessaire de toute façon si on utilise les données temps d'arrivée), modélisé par exemple par tracé de rais. En manipulant ces deux dernières propositions, il faut faire toutefois attention au fait qu'elles sont associées à un récepteur particulier et non pas au réseau de récepteurs (s'il existe) et à l'onde observée. Par conséquent, ces fdp a priori sont à prendre en compte au niveau des fdp-individuelles a posteriori : ap (X) et ati (X). Maintenant que les fdp-temps et fdp-polarisation sont normalisées correctement, on peut les combiner si nécessaire, ou bien travailler sur chacune individuellement, afin d'obtenir les incertitudes spatiales a posteriori liées à la localisation. Etant donné que nous avons les valeurs de ce (X) sur l'ensemble de la grille bornée, nous pouvons classer les noeuds par fdp-conjointe décroissante. Le premier noeud correspond à la meilleure localisation (déjà vu). Ensuite, en faisant une boucle sur les noeuds (classés), nous pouvons tout simplement sommer chacune des valeurs de la fdp-conjointe jusqu'à ce que la somme soit égale à certaines valeurs reconnaissables comme 0. 683, qui, dans le cas d'une fonction gaussienne, correspond à l'intégrale de cette fonction entre plus ou moins un écart type de sa valeur moyenne. L'ensemble des noeuds de la grille cible pour lesquels cette somme est inférieure à 0.683 appartiennent donc au domaine pour lequel il y a 68. 3% de chance que la localisation soit dans ce domaine. Il est alors possible d'analyser la distribution de l'ensemble de ces noeuds afin de quantifier l'incertitude sur la localisation par exemple à l'aide de la matrice de covariance, comme cela a été présenté auparavant. Nous avons illustré la méthode avec une valeur de 0.683, mais elle peut-être quelconque. La figure 5 est une représentation sous forme de courbes d'isoprobabilité du résultat d'un traitement conforme à celui qui vient d'être décrit. La fdp est définie de façon absolue sur l'ensemble des noeuds de la grille de la zone cible, ce qui permet d'obtenir des courbes d'isoprobabilité. Ces courbes permettent de visualiser l'incertitude spatiale a posteriori sur la localisation du séisme. APPLICATIONS DE L'INVENTION Comme on l'aura compris, de nombreuses applications de la technique proposée dans le domaine sismique sont envisageables. Par exemple celles qui nécessitent la determination de l'épicentre d'un (ou plusieurs) séisme (s), telles que la surveillance de sites présentant ou susceptibles de présenter de la simicité induite ou naturelle (mines, sites géothermiques, sites pétroliers, sites de stockage de gaz, etc. ) ; ou encore, les applications qui sont basées sur l'analyse de la localisation de séismes, par exemple, la cartographie de fractures actives ou des zones instables dans un réservoir pétrolier. Plus généralement, cette méthode peut être appliquée pour la localisation de tout autre type de source qui génère des ondes, ces ondes étant enregistrées sur un ou plusieurs récepteurs qui permettent d'en mesurer la polarisation. REVENDICATIONS : 1. Procédé de détermination d'une distribution de probabilités de localisation d'une source d'ondes observées par au moins un récepteur multi-axial, selon lequel : - on détermine à partir des signaux enregistrés sur les différentes composantes dudit récepteur la polarisation des ondes au niveau de celui-ci, - on détermine pour une pluralité de points d'une zone cible où la source est susceptible de se trouver une polarisation théorique attendue pour ces mêmes ondes audit récepteur pour le cas où la source se serait trouvée au point de la zone cible considérée, - on calcule pour les différents points considérés dans la zone cible une probabilité de localisation qui est fonction de la polarisation attendue au récepteur et de la polarisation mesurée. 2. Procédé selon la revendication 1, caractérisé en ce qu'on calcule pour les différents points considérés dans la zone cible une probabilité de localisation qui est fonction, d'une part, des composantes d'une matrice de covariance ou d'auto-corrélation calculée pour une onde donnée à partir des signaux enregistrés au niveau du récepteur et, d'autre part, des composantes du vecteur polarisation théorique attendu audit récepteur, pour cette même onde, et pour le cas où la source se serait trouvée au point de la zone cible considérée. 3. Procédé selon l'une des revendications 1 ou 2, caractérisé en ce que, la source d'ondes étant observée en une pluralité de récepteurs et/ou plusieurs polarisations d'ondes étant observées pour un même récepteur, on calcule pour les différents points de la zone cible considérés une probabilité de localisation qui est le produit des différentes probabilités de localisation pour les différents récepteurs et/ou les différentes polarisations d'ondes qui y sont observées. 4. Procédé selon l'une des revendications 1 à 3, caractérisé par le traitement suivant : - on détermine à partir des signaux relevés par un ou plusieurs récepteurs le ou les temps d'arrivée d'ondes au niveau de celui-ci, - on détermine pour les différents points considérés dans la zone cible un ou des temps d'arrivée théoriques attendus audit récepteur pour le cas où la source se serait trouvée au point de la zone cible considérée, - on détermine pour les différents points de la zone cible une probabilité conjointe de localisation qui est le produit d'une part d'une probabilité de localisation fonction des polarisations attendues et mesurées au (x) récepteur (s) et d'autre part d'une probabilité de localisation fonction des temps d'arrivée attendus et mesurés au (x) récepteur (s). 5. Procédé selon la revendication 2 prise seule ou en combinaison avec l'une des revendications 3 et 4, caractérisé en ce qu'une probabilité de localisation en un point X est égale à u, (X) = ( (X) exp- (X) (X) (X) 1 2 Â, Â, Px, Py et Pz : les composantes de la polarisation attendue exprimées dans la base propre de la polarisation observée, Â. i, Â. 2 et #3: les valeurs propres de la matrice de covariance calculée sur les données observées, i étant la plus grande de ces valeurs propres et À3 la plus petite, A étant un coefficient de normalisation dans l'espace, pp (X) étant une densité de probabilité a priori sur les données ou le modèle ; ou est un produit de cette probabilité avec d'autres probabilité. 6. Procédé selon l'une des revendications précédentes, caractérisé en ce qu'on détermine les noeuds de la zone cible pour lesquels la polarisation et/ou le temps mesurés sont compris entre plus ou moins un écart type autour de leurs valeurs modélisées, les autres noeuds étant considérés comme correspondant à une probabilité de localisation nulle. 7. Procédé selon la revendication 6, caractérisé en ce qu'on met en oeuvre une analyse statistique sur les points considérés dans la zone cible pour lesquels la probabilité de localisation n'est pas nulle. 8. Procédé selon la revendication 5 prise seule ou en combinaison avec l'une des revendications 6 ou 7, caractérisé en ce que la densité de probabilité a priori est : 9. Procédé selon la revendication 5 prise seule ou en combinaison avec l'une des revendications 6 ou 7, caractérisé en ce que la densité de probabilité a priori est fonction de la distance parcourue par l'onde considérée dans le milieu de propagation. 10. Procédé selon la revendication 5 prise seule ou en combinaison avec l'une des revendications 6 ou 7, caractérisé en ce que la densité de probabilité a priori est fonction du temps de propagation de l'onde considérée.