Brevet FR2864629A1 — publié le 20050701

Procédé de surveillance d'événements sismiques La présente invention concerne des procédés d'acquisition de données sismiques en surveillant passivement des événements sismiques tels que des signaux acoustiques générés par la production d'un réservoir d'hydrocarbure ou l'utilisation d'une stimulation hydraulique telle que la fracturation de couches de roche pour améliorer la production d'hydrocarbure d'un puits ou d'un réservoir. Plus spécifiquement, elle concerne de tels procédés utilisant des procédés sismiques pour déterminer les caractéristiques de source et l'emplacement de tels événements. CONTEXTE DE L'INVENTION La surveillance sismique est connue comme un procédé avec un horizon d'observation qui pénètre beaucoup plus profondément dans le réservoir d'hydrocarbure que n'importe quel autre procédé employé dans l'industrie pétrolifère. L'exploitation des procédés sismiques a été proposée en vue de la surveillance du réservoir. Dans la surveillance sismique conventionnelle, une source sismique, telle que des canons à air, des vibrateurs ou des explosifs, est activée et génère une énergie acoustique suffisante pour pénétrer dans la terre. Des parties réfléchies ou réfractées de cette énergie sont ensuite enregistrées par des récepteurs sismiques tels que des hydrophones ou des géophones. Dans la surveillance sismique passive, il n'y a pas de source activement commandée ou déclenchée. L'énergie sismique est générée à travers des événements dits micro-sismiques causés par des déplacements et changements souterrains qui entraîne au moins partiellement des ondes acoustiques qui peuvent à leur tour être enregistrées en utilisant les récepteurs connus. Sans parler du problème de détection des événements micro-sismiques souvent faibles, leur interprétation est difficile car ni l'emplacement de la source, ni la signature de la source, ni ses caractéristiques ne sont connus a priori. Cependant, la connaissance de ces paramètres est essentielle pour déduire d'autres paramètres de réservoir qui permettraient d'améliorer la commande du réservoir. Un domaine spécifique dans le secteur de la surveillance sismique passive est la surveillance de la fracturation hydraulique. Pour améliorer la production ou bien lorsque des réservoirs sont utilisés pour des raisons de stockage, les travailleurs de l'industrie du pétrole et du gaz exécutent une procédure connue sous le nom de fracturation hydraulique. Par exemple, dans des formations où le pétrole ou le gaz ne peuvent pas être facilement ou économiquement extraits de la terre, une opération de fracturation hydraulique est communément réalisée. Une telle opération de fracturation hydraulique comprend le pompage de grandes quantités de fluide pour induire des fissures dans la terre, créant ainsi des passages par lesquels le pétrole et le gaz peuvent circuler. Après qu'une fissure a été générée, du sable ou un autre matériau est habituellement ajouté à la fissure, de telle sorte que lorsque la terre se referme après que la pression soit relâchée, le sable aide à garder la terre écartée. Le sable fournit ensuite un passage conducteur pour que le pétrole et le gaz circulent à partir de la fracture récemment formée. Cependant, le processus de fracturation hydraulique ne fonctionne pas toujours très bien. Les raisons de cela sont relativement inconnues. De plus, les fractures hydrauliques ne peuvent pas facilement être observées, car elles sont typiquement situées à des milliers de pieds en dessous de la surface de la terre. Ainsi, les membres de l'industrie du pétrole et du gaz ont recherché des procédés de diagnostic servant à déterminer où les fractures se trouvent, de quelle dimension elles sont, jusqu'où elles vont et jusqu'à quelle hauteur elles se développent. Un appareil et un procédé de diagnostic destinés à mesurer la fracture hydraulique et la déformation de la roche autour de la fracture sont donc nécessaires. Dans des tentatives précédentes pour résoudre ce problème, certains procédés ont été développés afin d'établir des cartes des fractures. Par exemple, l'un de ces procédés implique la détection sismique. Dans une telle opération de détection sismique, des micro- tremblements de terre générés par la fracturation sont analysés par des appareils de mesure sismiques, par exemple, des accéléromètres. Une étude récente sur l'utilisation de l'imagerie micro-sismique pour la stimulation de fracture a été publiée par J.T. Rutledge et W.S. Phillips. Dans une configuration opérationnelle typique telle que décrite de manière plus détaillée à la Figure 1 ci-dessous, des géophones à trois composants ont été utilisés pour surveiller un puits pendant la fracturation. Les enregistrements des géophones ont été ensuite convertis en temps d'arrivée et emplacement de source en utilisant une méthode des moindres carrés itérative. La présente invention cherche à améliorer la quantité d'informations obtenues grâce à l'imagerie micro-sismique d'un réservoir dans le cas particulier des opérations de fracturation. RÉSUMÉ DE L'INVENTION L'invention décrit un procédé de traitement d'événements sismiques passifs comprenant des événements micro-sismiques ou la fracturation pour déterminer les caractéristiques de source, le temps origine ou l'emplacement de l'origine de ces événements au moyen d'une inversion de la forme d'onde. Contrairement aux procédés connus, le procédé de la présente invention peut être appliqué à la forme d'onde telle qu'elle est enregistrée et ne nécessite pas, par exemple de détection de phases sismiques spécifiques (telles que les ondes P ou S) ou d'autres paramètres dérivés des données (par exemple, les angles de polarisation). La forme d'onde complète est constituée de données enregistrées en utilisant des géophones à trois composants. De préférence, les signaux obtenus sont filtrés par filtre basse-bas ou par limitation de largeur de bande à une gamme de fréquence de 100 Hz ou moins, ou de préférence à 50 Hz ou moins. L'algorithme est adéquat pour l'inversion d'un milieu arbitraire hétérogène et profite d'un bon modèle de vitesse et de densité, s'il est disponible. Une version alternative de l'algorithme d'inversion (avec l'emplacement ou le temps origine de la source sismique déterminés indépendamment) peut être utilisée pour l'inversion, pour les caractéristiques ou le mécanisme de la source uniquement. Un exemple préféré d'une importante caractéristique de source est son tenseur de moment. L'algorithme utilise de préférence la réciprocité de la source et des récepteurs en évaluant des fonctions de Green dans un milieu arbitraire hétérogène depuis les emplacements de récepteur. Ces fonctions de Green sont ensuite inversées pour évaluer des sismogrammes synthétiques à cause d'un mécanisme de source arbitraire provenant des emplacements de source. En utilisant de préférence des algorithmes de recherche connus per se tels qu'une recherche de grille sur tous les emplacements de source et temps origines possibles, les sismogrammes synthétiques de forme d'onde complète sont adaptés aux donnés par la méthode des moindres carrés. L'estimation initiale du temps origine est définie par corrélation croisée de données et de synthétiques dus à un mécanisme de source arbitraire. Le temps origine inversé est déterminé par une recherche de grille autour de cette estimation initiale. L'algorithme est robuste par rapport au bruit blanc ajouté aux sismogrammes synthétiques et particulièrement adapté pour les données basse fréquence dans la bande de fréquence allant de 0 Hz à 100 Hz, de préférence de 0 Hz à 50 Hz. Ces aspects de l'invention ainsi que d'autres aspects sont décrits en détail dans les exemples suivants et les dessins qui les accompagnent. BRÈVE DESCRIPTION DES DESSINS L'invention va maintenant être décrite, à titre d'exemple uniquement, en référence aux dessins qui l'accompagnent, dans lesquels: La Figure 1 montre une illustration schématique d'une opération de fracturation; La Figure 2 est un organigramme d'étapes exécutées dans un exemple de la présente invention; et La Figure 3 est une comparaison de données synthétiques avec des données dérivées en utilisant un exemple de la présente invention. DESCRIPTION DÉTAILLÉE Une configuration opérationnelle typique destinée à surveiller une fracturation hydraulique est illustrée à la Figure 1 avec un puits de traitement 11 et des matrices de géophone 121, 131 situées dans des puits voisins ou des trous 12, 13. Pendant l'opération de fracturation, un fluide est pompé depuis la surface 10 dans le puits 11 ce qui fait se fracturer la formation alentour dans une couche de support d'hydrocarbure 101. Des ondes acoustiques 14 générées par la fracture 111 se propagent à travers la terre et sont enregistrées par les géophones à trois composants des deux matrices 121, 131. Pour la présente invention, il est supposé que trois composants de la valeur en fonction du temps de la vitesse acoustique des particules (ou déplacement des particules) en plusieurs (N_r) récepteurs de fond ont été enregistrés pendant une émission acoustique. De plus, il est supposé l'existence d'un modèle de vitesse (d'une complexité arbitraire) du volume de terre à travers lequel les ondes sismiques circulent. La qualité du modèle de vitesse peut être caractérisée par la longueur de l'intervalle de temps T i (i= 1. .N_r) pour lequel on estime qu'un sismogramme synthétique peut être adapté aux données. Ces intervalles de temps comprennent de préférence au moins l'arrivée d'onde S à tous les récepteurs. L'utilisation du déplacement des particules est préférée car il stabilise l'inversion alors que la vitesse acoustique des particules est plus oscillatoire que le déplacement des particules. Pour trouver les paramètres de source significatifs tels que le vecteur d'emplacement x_s, le temps origine t 0 et le tenseur de moment M, l'écart entre les sismogrammes synthétiques et les données est minimisé. Dans cette inversion, l'écart est défini par l'équation [1] : N,. 3 T, = E E (d(,4,t_0) U1(x,x,t,M))2dt i=0;=0 0 [1] Où d j indique un composant de la vitesse acoustique des particules enregistrée au 1eme récepteur et U_j est le jème composant du sismogramme synthétique au 1ème récepteur à cause d'une source située en x s caractérisée par un tenseur de moment M. Pour faciliter la description, les caractères suivant un trait de soulignement apparaissent en indice dans les équations. Les paramètres de source qui minimisent l'équation [1] comprennent la solution inversée. Le jeme composant d'un sismogramme synthétique au ieme récepteur xlr dû à des sources aux emplacements x s peut être évalué grâce à la relation bien connue suivante [2] (x,,xr,t,M.) = E Gki,rn (34 y Xr y t) * Mkm(Xs,t) Ici * est une convolution en temps, G kj,m est la dérivée de la fonction de Green le long du mème axe des coordonnées et M jk est un tenseur de moment d'une source point située en x_s. Le minimum des moindres carrés de l'écart donné par l'équation [1] est en général non unique. Pour alléger ce problème, on préfère émettre deux hypothèses: Tout d'abord, il y a approximation de la source comme une source point x s unique de telle sorte que la somme sur x s dans l'équation [2] disparaît. Deuxièmement, la fonction temps-source peut être approximée comme une fonction source delta de telle sorte que la convolution dans l'équation [2] est remplacée par une multiplication. En utilisant ces approximations l'équation [2] se réduit à uj (Xe, XÉr, t, m) = Gii,k (X8, x sr, t) M,k (Xa) Gi3,l (Xe, t) Xr, 'M11(x8) + G29,2(x8, X:. , t) M22(Xa) + G313(xe, xri, t) É M33 (x$) + (G2J,i (Xa, Xr, t) + G1?,2(x3, 7t., t)) É M21 (Xe) + (G3j,1(xi,7 t) +G3?,1(xe,Xrrt))) É M31(xi) + (G3121Xa, xr, t) + G3j,2(x3, X; , t)) É M32 (x5)' [3] Il est connu que l'équation [3] possède une seule solution pour M avec un temps origine t_0, un emplacement de source point x s et un modèle d'inversion fixes. Ainsi, les substitutions possibles parmi les paramètres de source peuvent être minimisées par une recherche de grille sur des emplacements sources et des temps origines pour obtenir les tenseurs de moments les mieux adaptés. La recherche de grille de tous les temps origines possibles est numériquement coûteuse et est donc accélérée en estimant le temps origine par une corrélation croisée des synthétiques et des données puis en utilisant la recherche de grille autour de l'estimation initiale. Le procédé utilisé comprend les étapes suivantes, telles que l'illustre la Figure 2: - Suivre un enregistrement de données acoustiques à partir d'une fracture (Étape 20) ; - Estimer le temps origine initiale t 0(x s) à chaque emplacement de source possible (x s) (Étape 21) ; - Exécuter une recherche de grille autour du temps origine estimée pour chaque emplacement de source (Étape 22). Pour chaque temps origine, trouver la solution unique M(x s,t 0(x s)) (minimum des moindres carrés) (Étape 23) et évaluer l'écart des moindres carrés entre les données et les synthétiques (Étape 24) ; et - Stocker la solution la mieux adaptée pour chaque 15 emplacement de source (Étape 25). Le tenseur de moment de fracture avec le temps origine et l'emplacement peuvent être également évalués (Étape 26) comme cela est décrit cidessous pour trouver les caractéristiques de la fracture. L'estimation initiale du temps origine est évaluée par corrélation croisée des données et des sismogrammes synthétiques pour un mécanisme de source choisi, par exemple un rejet vertical. La corrélation croisée est évaluée sur l'intervalle de temps $(0, T _j) pour chaque récepteur j. Les valeurs absolues des composants correspondants pour chaque récepteur sont soumises à une corrélation croisée et les décalages de temps de la corrélation croisée maximale pour chaque composant sont calculés. L'utilisation des valeurs absolues des sismogrammes pour la corrélation croisée réduit la dépendance du mécanisme de source inconnu. Les décalages de temps de chaque composant et les temps origines connus des sismogrammes synthétiques permettent une estimation du temps origine absolu t ij pour chaque composant i et récepteur j. Les estimations sont pondérées par l'amplitude maximum des sismogrammes enregistrés pour réduire les mauvaises estimations résultant des traces de corrélation croisée dominées par le bruit. Il vaut la peine de remarquer que l'utilisation de l'amplitude maximum en guise de poids dans le calcul de la moyenne du temps origine suppose que le rapport signal-bruit est proportionnel à l'amplitude maximum des sismogrammes enregistrés. L'estimation finale du temps origine est donc une moyenne pondérée arithmétique avec les poids des amplitudes maximum A ij du i ème composant au ème composant: [4] r0E=0ta t0(x) - EN3 Cette corrélation croisée peut être encore améliorée aux dépens d'un calcul plus intensif dans le temps en utilisant des enveloppes de signaux à la place des amplitudes. Le véritable temps origine est ensuite trouvé par recherche de grille autour de l'estimation initiale du temps origine à l'intérieur de la période dominante [la plus courte] dans le signal. La limitation de la recherche de grille à la période dominante du signal nécessite que l'estimation initiale du temps origine [4] soit située à l'intérieur de la période dominante. C'est typiquement le cas pour l'arrivée de l'onde S. La recherche de grille autour de l'estimation initiale du temps origine [4] élimine les problèmes liés aux sauts de cycle car la fonction de corrélation croisée tend à atteindre son apogée toutes les demi-périodes de la période dominante (habituellement, la période minimum présente dans les données). La longueur de l'échelon de temps dans la recherche de grille est définie de manière à obtenir la précision souhaitée de l'écart [1] . En supposant que les sismogrammes synthétiques correspondent aux données (c'est à dire en utilisant le mécanisme de vrai moment et en évaluant les sismogrammes synthétiques dans le vrai modèle à partir du vrai emplacement de source), un écart normalisé d'un signal harmonique avec une période T, dû à un décalage de temps aT dans le temps origine, peut être évalué de la manière suivante [5] fT[$jfl(wt) sin (w t + c 7) )]2dt =1 _ E 2 fTo [sin (cd; t)]2 t La définition d'erreur dans l'équation [5] a un maximum de 2 pour 1/2 décalage de période et même un décalage de temps cause une grande erreur pour un écart défini de manière analogue à l'équation [1]. La longueur de l'échelon de temps pour la recherche de grille peut être définie à 2aT, pour lequel l'erreur maximum d'évaluation de l'écart atteint une certaine limite. Par exemple, un décalage de 0, 05T$ (a=0,05) peut causer une erreur relative E=0,05. Ainsi, une recherche de temps origine avec un échelon de grille de O,1T (T est la période dominante dans les sismogrammes) ne devrait pas causer une erreur d'évaluation dans la fonction d'écart supérieure à 0,05. La dernière partie du procédé consiste à identifier une solution unique M(x s, t 0(x s)) pour chaque temps origine et emplacement de source. Il est connu que le tenseur de moment avec l'ajustement minimum des moindres carrés de l'équation [1] est: [6] (x8) (A iii(x.5)Di(xs). Ici, M I(barre) est le lem" composant de six éléments vecteurs: M I(barre) = M 11, M 2(barre) = M 12 = M 21, M 3(barre) = M 22, M 4(barre) M 13 = M 31, M 5(barre) _ M 3 = M32, M 6(barre) = M 33, et D possède six éléments indépendants [7] N,. 3 T; E E Dk (xs) gik(xs, xf, t -- to)d1(xi., t)dt. i=O j =o Ici, k=0...5 et g jk suivante: [8] est défini par la notation 911(x8, xr t) 912(xs,xi.,t) 913{x8:xr, t) 914(x8, xi., t) 9j5(xs, xr, t) 916(xs, xr, t) G11,1(Xie, xr, t) G21,1(x8 xr, t) + G11,2(xs 7 xr, t) G21, 2 (x8 2 xr t) G31,1(xs, xr, t) + G11,3{x8, xr, t) G31s2(xs, xr, t) + G21, 3(x8, x, , t) G31,3(xs, xr, t) Enfin, A est une matrice 6 x 6 avec les éléments: [9] 13 N,. 3 f T, Akj (x$) = E E g (x3, xT, t)9a (x$, x;,, t) dt. i=0.1-o Les étapes d'intégration de [7] et [9] peuvent être accélérées en utilisant une fenêtre de temps t min à t max, où t min est un temps d'arrivée d'une première énergie provenant de la source (fracture) telle qu'identifiée par un détecteur d'événement et t max est le temps maximum pour lequel les formes d'onde sont associées, par exemple le temps d'arrivée de la phase avec l'amplitude maximum. Cette modification exclut l'effet de réflexions ou d'ondes de tube dans les données enregistrées. Lors de l'extraction du tenseur de moment M à partir des enregistrements à trois composants du champ d'ondes en résolvant le problème d'inversion des moindres carrés, il se peut que la solution ne soit pas stable car par exemple la matrice A peut avoir un problème de rang. Pour obtenir une solution stable de ce problème, une régularisation algébrique peut être appliquée. Pour régulariser le problème, seules les valeurs propres les plus grandes sont sélectionnées avec un numéro de conditionnement inférieur à une limite prédéfinie et une décomposition tronquée des valeurs singulières est réalisée. Le degré de singularité de la matrice est mesuré en calculant le numéro de conditionnement de la matrice pour chacune des valeurs propres. Le numéro de conditionnement est exprimé par le rapport entre chaque valeur propre et la valeur propre la plus grande. Le critère de seuil consiste à vérifier que le numéro de conditionnement ne dépasse pas la valeur de seuil. Chaque numéro de conditionnement est comparé à la valeur de seuil. Le nombre de valeurs propres qui satisfont au critère de seuil est équivalent au rang de la matrice. Une fois que le nombre de valeurs propres k qui fournissent des solutions indépendantes linéaires est déterminé, une décomposition de valeur singulière tronquée est utilisée pour résoudre le problème inverse. La nouvelle solution inverse est calculée par l'expression suivante: [10] k E ul D ' vi i=l ai Où M barre est le tenseur de moment stabilisé, D est le vecteur de données, u et v sont les vecteurs propres et o sont les valeurs propres obtenues par la décomposition de valeur singulière. Dans l'équation [10] seuls les vecteurs propres correspondant aux valeurs propres k acceptables sont utilisés pour inverser la matrice. Il est également faisable d'associer avec chaque dispositif d'enregistrement ou trace une fonction de pondération qui indique la qualité du récepteur et/ou des données enregistrées. Ces poids pourraient être introduits dans les présentes équations [7] et [9]. La fonction de Green synthétique dans l'équation [3] est ensuite évaluée en calculant trois fois des simulations de forme d'onde complète N r (en utilisant une différence finie). Pour chaque récepteur à trois composants, trois réponses dues à trois sources de force orthogonale simple aux positions des récepteurs sont calculées et des dérivées de la vitesse (ou du déplacement) sont stockées à chaque emplacement de source Mk possible, x s. Les sismogrammes synthétiques sont évalués avec une fonction delta comme fonction source. En utilisant la réciprocité, des dérivées des fonctions de Green pour chaque emplacement de source possible à chaque position de récepteur sont évaluées. L'équation [3] montre que six traces à chaque emplacement possible de source doivent être stockées. L'équation ci-dessus fournit une série complète d'étapes pour calculer le tenseur de moment M d'enregistrements à trois composants du champ d'ondes. Le tenseur lui-même est alors décomposé pour fournir des paramètres caractéristiques de la fracture. Des procédés pour décomposer le tenseur de moment M ont été développés dans le but d'analyser des tremblements de terre et sont décrits par exemple par V. Vavrycuk dans: Journal of Geophysical Research, Vol 106, N B8, août 10, 2001, 16,339-16,355. Les paramètres obtenus par une telle décomposition comprennent les paramètres normaux de la fracture n, la direction de glissement N, et des produits des coefficients de Lame avec le glissement u de la fracture, à savoir, pu et Àu respectivement. Alternativement, le tenseur de moment peut être inversé pour un groupe de paramètres comprenant l'orientation des axes de pression P et de tension T, le paramètre K = À/u et l'inclinaison a du glissement u depuis la fracture. Ces paramètres fournissent des informations sur l'orientation de la fracture et la direction de glissement qui peuvent à leur tour être utilisées pour contrôler l'opération de fracturation hydraulique. La précision de l'inversion des données enregistrées d j au tenseur de moment M de la source peut être encore améliorée en limitant la largeur de bande de la fréquence des données. Tandis que le fai-: de limiter les données à une gamme de fréquence comprise entre 0 et 100 Hz permet l'obtention de résultats satisfaisants, une précision améliorée est gagnée en limitant encore plus les données à une gamme de fréquence comprise entre Oet 75Hz et même à une gamme de fréquence comprise entre 0 et 50 Hz. À la Figure 3, il est montré un relevé de mesures de vitesse 31 d'un géophone (synthétique) dans les directions x, y et z superposées avec les traces correspondantes 32 ]0 recalculées en utilisant le tenseur de moment dérivé par le procédé décrit ci-dessus (avec un modèle de vitesse connu). Le procédé décrit ci-dessus et les variantes de ce dernier peuvent être appliqués à l'analyse de n'importe 15 quel autre événement micro-sismique. REVENDICATIONS 1. Procédé de surveillance passive d'un emplacement souterrain comprenant les étapes d'obtention de signaux multi-composants d'un événement microsismique à l'intérieur de l'emplacement; et de réalisation d'une inversion de forme d'onde pour déterminer des paramètres représentant les caractéristiques de source dudit événement micro-sismique. 2. Procédé selon la revendication 1 dans lequel les enregistrements de signaux sont, au moins dans le but de déterminer les caractéristiques de source, filtrés par filtre passe-bas ou limitation de largeur de bande à une gamme de fréquence située entre 0 et 100 Hz. 3. Procédé selon la revendication 1 dans lequel l'événement microsismique est causé par une opération de fracturation dans un puits de forage. 4. Procédé selon la revendication 1 comprenant l'étape d'évaluation d'une fonction de Green pour dériver les caractéristiques de source à partir des signaux obtenus. 5. Procédé selon la revendication 1 dans lequel les signaux obtenus sont traités pour identifier les événements d'onde P et d'onde S avant l'inversion du champ d'ondes. 6. Procédé selon la revendication 1 dans lequel les paramètres des caractéristiques de source sont représentés par un tenseur de moment et/ou un emplacement de source et/ou un temps origine. 7. Procédé selon la revendication 1 comprenant également l'étape d'utilisation de décomposition de valeur unique pour stabiliser l'inversion de forme d'onde. 8. Procédé selon la revendication 1 comprenant l'étape de minimisation de la différence entre les signaux obtenus et les signaux synthétiques. 9. Procédé selon la revendication 1 comprenant l'étape de minimisation de la différence entre les signaux obtenus et les signaux synthétiques avec les signaux synthétiques dépendant des caractéristiques de source estimées. 10. Procédé selon la revendication 9 dans lequel l'étape de minimisation de la différence entre les signaux obtenus et les signaux synthétiques comprend une recherche sur les emplacements de source et les temps origines pour trouver des caractéristiques de source estimées. 11. Procédé selon la revendication 9 dans lequel l'étape de minimisation de la différence entre les signaux obtenus et les signaux synthétiques comprend une recherche de grille sur les emplacements de source et les temps origines pour trouver des caractéristiques de source estimées. 12. Procédé selon la revendication 10 comprenant en outre les étapes de estimation du temps origine initial en des emplacements de source possibles; exécution d'une recherche autour du temps origine estimé pour chaque emplacement de source; - pour lesdits temps origines, découverte de l'unique solution d'un tenseur de moment de la source; - évaluation de l'écart mesuré par fonction des moindres carrés entre les signaux enregistrés et les signaux synthétiques dérivés par calcul des signaux causés par une source dudit tenseur de moment aux emplacements du récepteur; et ]0 - stocker la solution la mieux adaptée pour chaque emplacement de source. 13. Procédé selon la revendication 1 dans laquelle une fonction source de la fracture est approximée par une 15 fonction delta.