Brevet FR2958411A1 — publié le 20111007

Procédé d'analyse spectrométrique et dispositif apparenté L'invention concerne un procédé d'analyse spectrométrique et un dispositif apparenté. L'invention s'applique notamment aux domaines de la spectrométrie Gamma, de la spectrométrie X et de la spectrométrie de masse à temps de vol. L'analyse de données spectrométriques a notamment pour objectif l'identification d'un nombre inconnu de pics superposés à un fond également 1 o inconnu, lesdits pics étant représentatifs de quantité physiques présentent dans le milieu à analyser. Plusieurs difficultés affectent cette analyse. La forme des pics présents dans le spectre n'est habituellement modélisable que partiellement par des fonctions paramétrables. De la même façon, il est difficile de 15 décomposer les fonds en utilisant une base de fonctions paramétrables. Le fait que le nombre de pics soit inconnu complique également l'analyse. Plusieurs méthodes d'analyse automatique de spectres sont proposées dans l'état de la technique pour contourner ces difficultés et faciliter l'analyse de données spectrométriques. Deux catégories des 20 méthodes peuvent être distinguées. Une première catégorie comprend les méthodes d'analyse s'appuyant sur des bases de données. Ces techniques sont répandues dans le domaine de la spectrométrie nucléaire et de la spectrométrie de masse. Le 25 principe est de rechercher les pics au sein d'une base de données qui les caractérise en position et en intensité. Ce type de méthode permet de chercher les positions et intensités sur des espaces continus en fixant le nombre de composantes. Il est alors possible d'estimer la quantité de chaque élément composant la base de donnée en minimisant une fonction de coût. 30 La forme des pics est dans ce cas supposée déterministe et connue. Pour ce qui est du fond, celui-ci est modélisé au préalable sous chaque pic d'une région d'intérêt donnée par une fonction polynomiale. L'utilisation d'une procédure d'ajustement et de modèles paramétriques figés sur des données bruitées induit une erreur de reconstruction constituée de deux composantes, 35 la première correspondant à une erreur statistique et la seconde à une erreur de modèle. De plus, cette approche suppose l'exactitude de la base de données de référence, cette hypothèse étant rarement vérifiée dans la mesure où une part expérimentale importante contribue à l'élaboration desdites bases. Il est alors difficile d'estimer l'incertitude totale liée à l'analyse. Une seconde catégorie comprend des méthodes d'analyse basées sur une détection automatique de pics. Cette approche ne requiert pas l'utilisation de bases de données mais suit une approche statistique fondée sur la minimisation d'un risque ou la maximisation de vraisemblance. 1 o La recherche de la position desdits pics est effectuée sur un espace discret constitué des canaux du spectre acquis. En effet, des analyseurs multicanaux sont habituellement utilisés pour analyser des spectres tels que des spectres X ou y. Des contraintes techniques, comme la résolution limitée des convertisseurs analogiques-numériques par exemple, 15 ne permettent pas de générer directement un spectre continu. C'est pourquoi ces analyseurs présentent en sortie des histogrammes composés de plusieurs canaux. Un canal j d'histogramme correspond par exemple à un intervalle d'énergie B1 , les quantités à mesurer, des photons par exemple, sont classés par canaux en fonction de leur niveau d'énergie. Deux pics 20 distincts peuvent être détectés si ils sont séparés d'au moins un canal de spectre. Ces canaux sont également appelés canaux d'accumulation dans la suite de la description. L'ensemble des canaux d'accumulation représentatifs d'un spectre composent un histogramme. Dans le cadre des méthodes d'analyse basées sur une détection 25 automatique de pics, l'estimation de l'intensité des pics s'opère sur un espace continu. Le nombre de pics peut être supposé connu à l'avance sur une région d'intérêt limitée ou bien estimé à partir des données. Les fonds peuvent être estimés empiriquement ou modélisés suivant des combinaisons linéaires de splines telles que les B-splines, par exemple. Dans tous les cas, 30 le nombre de pics est fini. Le nombre fini de pics recherchés et la forme paramétrée desdits pics induit une erreur statistique et une erreur de modèle lors de la reconstruction du spectre. Moyennant une hypothèse de bruit, ledit bruit correspondant à la nature statistique du bruit d'un canal d'histogramme, il est possible d'accéder à l'incertitude des mesures en utilisant l'inverse de la 35 matrice du Hessien. Cette estimation est efficace pour une hypothèse de bruit gaussien, mais n'est pas représentative de la réalité en cas de bruit Poisson, et ne prend pas en compte les éventuelles erreurs de modèle. Un exemple de méthode de détection automatique de pics est décrit dans l'article de S. Gulam Razul, W. J. Fitzgerald et C. Andrieu intitullé Bayesian mode/ select/on and parameter estimation of nuclear emission spectra using RJMCMC, Elsevier, Nuclear Instrument and Methods in Physics Research, A 497, pages 492-510, 2003. Cette approche est paramétrique au sens où un nombre de pics inconnu mais fini est considéré. Les méthodes d'analyse spectrométrique décrites précédemment ne 1 o permettent pas une détection fiable des pics présents dans le spectre lorsque le nombre de pics présents est potentiellement infini. De plus, il n'est pas possible d'estimer le niveau d'incertitude de ladite détection par canal d'accumulation. 15 Avantageusement, la méthode proposée prend en compte le bruit de quantification du au fait que les systèmes d'acquisition traditionnels mémorisent les échantillons de mesures dans des histogrammes. L'invention a aussi comme avantage de permettre une implémentation ne requérant pas un puissance de calcul important, les traitements étant faits par canaux 20 d'accumulation et non par échantillon. Un but de l'invention est notamment de pallier les inconvénients précités. A cet effet l'invention a pour objet un procédé d'analyse de 25 mesures spectrométriques, lesdites mesures C = {c1,...,cY} étant classées dans un histogramme, ledit histogramme étant composé de canaux d'accumulation, un canal j correspondant à un intervalle d'énergie Bi . Ledit procédé comporte au moins trois étapes de traitement. Une première étape détermine par canal d'accumulation j des distributions p représentatives de 30 l'amplitude et Z représentatives de la position des impulsions de Dirac composant le spectre normalisé de pics ainsi que des distributions q représentatives de l'arbre de Pôlya caractérisant le spectre normalisé de fond, lesdits éléments étant obtenus par échantillonnage de Gibbs. Une seconde étape de détecte des canaux significatifs, ladite détection étant 35 réalisée par application d'un test de Kolmogorov-Smirnov par canal d'accumulation de l'histogramme sur la base des résultats de la première étape. Une troisième étape identifie des régions significatives du spectre par regroupement des intervalles B. des canaux significatifs identifiés lors de la seconde étape, lesdites régions comportant au moins un pic significatif. Procédé selon un aspect de l'invention, le spectre de pics est déterminé pour tout canal d'accumulation d'indice j, en utilisant des estimateurs Monte-Carlo définis par une expression telle que : n L E(wnII1 IC) =-1 L dans laquelle : le correspond à la masse de probabilité du canal B . pour le 1 ième tirage du processus MCMC ; L correspond au nombre de tirages Monte Carlo de l'échantillonneur de Gibbs ; we [0,1] est la probabilité du spectre de pics dans le spectre complet défini comme étant la somme pondérée du spectre normalisé de pics et du spectre normalisé de fond ; w(l) représente le poids du spectre de pics pour le l ième tirage du processus MCMC ; n est le nombre de photons enregistré. Le spectre de fond est déterminé, par exemple, pour tout canal d'accumulation d'indice j , en utilisant l'expression suivante : E((1ùw)n(D I C) = n 1(1-L ~-~ dans laquelle : cl)(» représente les valeurs générées de la variable aléatoire cl)j à l'itération l de l'échantillonneur, J représentant la masse de probabilité du canal B. du spectre normalisé de fond pour chaque j L'incertitude associée à chaque canal énergétique est estimée par exemple en définissant un intervalle de crédibilité à 95% [1_,ri+] déterminé en utilisant les expressions suivantes : = min(Pr(w n n i _ 0,025) - n{w~~' Il i }T (L0,025L]) ~7+ = min(Pr(w n II1 ù 0,975) n{w»ri y }T (L0,975L]) dans lesquelles : pour tout x , la notation {x}T représente la collection des échantillons x triée par ordre croissant ; pour tout x et pour tout j la notation {x}T (j) représente le j ième échantillon de la collection x triée par ordre croissant. ainsi, {w(l)fl »}T représente la collection des échantillons générés w("IÉ triés par ordre croissant ; L.] indique la partie entière ; n représente le nombre total de photons enregistrés ; 17 représente une variable caractérisant l'intensité du spectre de pic pour le canal d'énergie considéré. Une quantité DL(F'J,EJ) est déterminée par exemple pour application du test de Kolmogorov-Smirnov (210) en utilisant l'expression : DL(F'J,E. ) = Lsup 2 n CDF ~(1 ~} (r7) ù CDF dans laquelle : E; correspond à la probabilité du canal Bi dans le spectre complet normalisé ; F. est une quantité correspondant au fond seul ; {E(;) } représente la collection des échantillons générés 1:1) pour lL; {Fr} représente la collection des échantillons générés Fi» pour lL ; CDF (I)l(r7)la fonction cumulative empirique de la distribution de Ei telle que CDF EY)} (17) = 1 (E») 17 représente une variable caractérisant l'intensité du spectre complet pour le canal d'énergie considéré. La présence d'un élément de pic significatif sur un intervalle B. est détectée par exemple si DL(FJ,EJ)> Ka, Ka étant défini tel que Pr(K dans lesquelles : la fonction 1(Cond) =1 si la condition Cond est vraie et 1(Cond) = 0 sinon ; {R, n} la collection de régions significatives du spectre ; Zk est la position de la composante k du processus de Dirichlet produit par la procédure MCMC ; Le procédé selon l'invention comporte par exemple une étape de 20 calcul d'un intervalle de crédibilité à T95% pour 2m en utilisant l'expression : C195% (2m) [{2 }T (L0,025L]),{2 }T (LO,9 / 5L])] dans laquelle la collection {2m' }T correspond à la collection d'échantillons {2.,`' } ordonnée par ordre croissant. 25 Le procédé selon l'invention comporte par exemple une étape d'estimation du centroïde ,ûmde chaque région de pics significatifs en utilisant l'expression suivante : L Diu) P. L I1(Nmi) > o) c=1 1=1 30 dans laquelle : ,un.,'» est une quantité estimée sur chaque région significative Rm et pour chaque tirage de la procédure MCMC en utilisant l'expression : \k=1 7 = I pk ("eI(Z" E Rm ) k =l le dénominateur correspond aux nombre de tirages MCMC pour lesquels il y a au moins une composante appartenant à la région Rm , les autres tirages ne pouvant pas être pris en compte dans l'estimation Monte Carlo. Le procédé selon l'invention comporte par exemple une étape de 1 o calcul d'un intervalle de crédibilité à 95% pour ,um en utilisant l'expression : CI95% (J` m) = L{ }T (L0,O25LJ),{,um`) }T (Lo,975L])] dans laquelle la collection {,um`) }T correspond à la collection 15 d'échantillons {,u, n(»} ordonnée par ordre croissant. Le procédé selon l'invention comporte par exemple une étape d'estimation 6m de l'écart-type am de la distribution d'énergie restreinte à la région Rm en utilisant l'expression : L 1=1 = 6m L 1I(Nm~) > 0) 1=1 dans laquelle : 20 N (l) 1(Z (l))2 1(z (l) = k k `\ k \k=1 N Rm) _ I p(1)Z(1)I(z ) 1 25 Le procédé selon l'invention comporte par exemple une étape de calcul d'un intervalle de crédibilité à 95% pour am en utilisant l'expression : CI95% ( m ) = [{o-,, }î (LO,025L]),{ }T (Lo,975L])J dans laquelle la collection {0 }T correspond à la collection d'échantillons {a,,( } ordonnée par ordre croissant. Le procédé selon l'invention comporte par exemple une étape d'estimation de la largeur de la région déconvoluée en sa base, définie 5 comme l'intervalle de plus forte probabilité à 99% de la restriction à Rm de la densité d'énergie (notée Am,99%) , en utilisant l'expression 1 o dans laquelle : A(m),99% = [QR12 (0,005),QR,) (0,995)] où QRi)(a)= CDF(Ri?~ )= a et CDF" ( ) = l Pi l) I~Zki~ 20 dans laquelle la collection {o(m",99%}1 correspond à la collection d'échantillons {O(m'99%} ordonnée par ordre croissant. L'invention a aussi pour objet un dispositif d'analyse de mesures spectrométriques, lesdites mesures C = {c1,...,cY} étant classées dans un 25 histogramme, ledit histogramme étant composé de canaux d'accumulation, un canal j correspondant à un intervalle d'énergie B1 , ledit dispositif étant caractérisé en ce qu'il comporte des moyens pour mettre en oeuvre le procédé décrit précédemment. D'autres caractéristiques et avantages de l'invention apparaîtront à l'aide de la description qui suit donnée à titre illustratif et non limitatif, faite en regard des dessins annexés parmi lesquels : - la figure 1 illustre le principe de la méthode d'analyse de spectre selon l'invention ; - la figure 2 donne un exemple de méthode de construction et d'analyse de spectre échantillonné ; - la figure 3 donne un exemple de mise en oeuvre du test de Kolmogorov-Smirnov dans le cadre d'une application de spectrométrie Bayésienne semi paramétrique ; - la figure 4 donne un exemple de traitement de construction séquentielle des régions de pics significatifs ; - la figure 5 illustre les analyses pouvant être conduites après application du test Kolmogorov-Smirnov dans le but d'analyser les pics détectés ; La figure 1 illustre le principe de la méthode d'analyse de spectre selon l'invention. Celui-ci comprend un module échantillonneur de Gibbs 100 ainsi qu'un module de construction et d'analyse de spectre échantillonné 101. Dans l'article de E. Barat et al. T. Dautremer et T. Montagu intitulé Nonparametric bayesian inference in nuclear spectrometry, IEEE Nuclear Science Symposium Record, NSS/Ml 2007, une extension des méthodes de détection automatique de pics décrites précédemment est présentée. Cette extension propose de considérer le nombre de pics comme potentiellement infini. Le principe est de rechercher une mesure de probabilité qui, par convolution d'un noyau donne une distribution de probabilité dont des échantillons rangés forment l'histogramme observé, un noyau étant ici une fonction positive, normalisée à 1, correspondant à la distribution du bruit instrumental variant avec la position sur l'histogramme, ladite position correspondant par exemple à un niveau d'énergie ou de masse. La nature statistique du bruit d'un canal d'accumulation est prise en compte car l'histogramme correspond au comptage rangé d'une collection de réalisations générée par une distribution de probabilité. Cette approche présente un caractère non paramétrique dans la mesure où une quantité pouvant évoluer dans tout l'espace des distributions de probabilités sur l'espace des positions est estimée. Ainsi, une éventuelle imperfection de modèle dans la forme des pics, correspondant à du bruit instrumental, sera prise en charge par la mesure de probabilité recherchée au moyen de la présence potentielle d'une infinité de pics. Par ailleurs, le fond est également modélisé par une mesure de probabilité obtenue par convolution avec le même noyau. La discrimination entre la mesure de probabilité correspondant aux pics et celle correspondant 1 o au fond est basée sur le fait que la distribution des pics sans convolution présente un caractère discret par rapport à la distribution sans convolution du fond, supposée continue à ce niveau de modélisation. D'autre part, cette méthode présente également un caractère paramétrique car la densité de la distribution de bruit instrumental est supposée paramétrique, mais ne 15 cherche pas nécessairement à estimer complètement la forme potentiellement complexe du pic. Cette densité peut être, par exemples, gaussienne pour un spectre gamma, ou correspondre à une fonction de Voigt pour un spectre X. Pour ces raisons, cette approche peut être qualifiée de semi-paramétrique. Cette méthode, à la différence des techniques 20 rappelées précédemment, garantit que la totalité des coups présents dans le spectre observé se retrouvent dans la somme du spectre et du fond déconvolué. Il y a ainsi conservation de l'intégrale du spectre. En partant d'un histogramme, un spectre déconvolué est obtenu avec d'un côté les pics et de l'autre le fond. L'intégrale dudit spectre est rigoureusement identique à celle 25 de l'histogramme de départ. La résolution du spectre obtenu est significativement améliorée par rapport au spectre initial, sans pour autant avoir limité la recherche à un nombre fini de pics Ce résultat est obtenu en utilisant la méthode d'échantillonnage de Gibbs, ladite méthode faisant partie des méthodes MCMC, acronyme venant 30 de l'expression anglo-saxonne « Markov Chain Monte Carlo ». La méthode proposée dans le cadre de l'invention utilisa un échantillonneur de Gibbs 100, ledit échantillonneur ayant notamment trois objectifs. Le premier objectif est de séparer automatiquement les pics du 35 fond sur l'ensemble du spectre, et ce de manière non-paramétrique. Le second objectif est d'améliorer la résolution en sortie en réalisant une déconvolution de forme gaussienne. En effet, les pics à détecter sont habituellement étalés avant traitement et ont une forme pouvant s'apparenter à une courbe gaussienne. La déconvolution permet de réduire la largeur des pics à détecter. Par conséquent, la largeur des pics est réduite et la précision de détection est ainsi améliorée. Enfin, l'échantillonneur fournit un tirage aléatoire des éléments constitutifs du spectre normalisé de pic et du spectre normalisé (à 1) de fond. Ces différents tirages aléatoires correspondent à : - l'amplitude p des impulsions de Dirac, lesdites impulsions étant également appelées masses de Dirac ; - la position Z desdites masses de Dirac ; - l'amplitude q du fond. p et Z sont des tirages aléatoires représentatifs du spectre normalisé des pics et sont exprimés dans l'espace des énergies sous forme de vecteurs. q est un tirage aléatoire représentatif du spectre normalisé de fond et exprimé sous forme de vecteur. L'échantillonnage de Gibbs appliqué dans le cadre de la spectrométrie tel que décrit dans l'article d'E. Barat et al. cité précédemment peut être mis en oeuvre afin d'obtenir les distributions de p, Z et q. Le module d'échantillonnage 100 est capable de traiter des échantillons de mesures collectés sous forme d'histogrammes 104 par un analyseur multicanaux, lesdites données d'entrée pouvant s'écrire sous la forme d'un vecteur C : C = {c~ , ... , cY } dans lequel les éléments ci représentent le nombre de coups dans le 1 ème canal d'accumulation B ., un nombre de coups c. étant équivalent à un nombre d'impulsions générées par le détecteur. D'autres données d'entrée sont requises par l'échantillonneur de Gibbs afin de générer les tirages aléatoires p, Z et q. Ainsi, des a priori, c'est-à-dire des distributions de probabilité conditionnelles, sont calculés en utilisant des techniques non paramétriques d'inférence Bayésienne. Le processus de Dirichlet et des arbres de Pôlya peuvent ainsi être utilisés pour déterminer lesdits à priori tels que décrits dans l'article de E. Barat et al. cité précédemment dans la description. La technique de déconvolution et de séparation des pics et des fonds présentée précédemment présente plusieurs propriétés. 1 o Une première propriété est que cette technique est non-paramétrique. En d'autres termes, la solution est recherchée dans l'ensemble contenant toutes les distributions de probabilités à valeurs réelles positives. Cet ensemble est noté M(9Z+). Une seconde propriété est liée à son caractère Bayésien. La 15 distribution a posteriori des processus de Dirichlet est disponible après échantillonnage de Gibbs et peut être déterminée par la collection de tirages aléatoires de p et Z . La distribution a posteriori des arbres de Pôlya peut être déterminée par la collection de tirages aléatoires de q. Le caractère non paramétrique de cette technique permet une 20 grande flexibilité de modélisation puisqu'elle s'affranchit de toute liste finie et déterminée à l'avance de positions possibles de pics d'énergie. Avantageusement, chaque tirage aléatoire intervenant dans la procédure d'échantillonnage MCMC 100 d'un spectre de pics déconvolués (processus de Dirichlet) est ainsi une distribution de probabilité discrète dont le nombre 25 de masses, c'est-à-dire de pics déconvolués, est potentiellement infini. Cette flexibilité de modélisation induit cependant une difficulté d'interprétation. Au cours des itérations de la procédure MCMC, il est possible que des masses de Dirac apparaissent en toute valeur d'énergie, conduisant asymptotiquement à une quantité mesurée quasi-continue. Afin 30 de fournir à l'utilisateur une liste de positions et d'intensités de pics pouvant être interprétée, il est possible de déterminer des régions significatives, c'est-à-dire des intervalles en énergie non explicable statistiquement par une fluctuation aléatoire du fond. Le caractère Bayésien offre ici un élément clé. L'accès aux distributions a posteriori permet en premier lieu de 35 mettre en oeuvre un estimateur du spectre de pics déconvolué mais permet également d'estimer l'incertitude associée à ce spectre de pics pour tout intervalle (région) en énergie. L'incertitude peut également être déterminée pour le spectre de fond. Dans la suite de la description, la génération des tirages aléatoires 5 par l'échantillonneur de Gibbs est aussi désignée par l'expression « procédure MCMC ». Ainsi, suite à la génération des tirages aléatoires par la procédure MCMC, l'objectif poursuivi est d'analyser un spectre échantillonné en canaux d'énergies dans le but de détecter et de caractériser des pics 1 o significatifs, c'est-à-dire des pics du spectre représentatifs de quantités physiques présentes dans le milieu mesuré par spectrométrie. Ainsi une détection des canaux significatifs est conduite 101, les canaux significatifs étant des canaux comportant des mesures relatives à un pic significatif. Les canaux significatifs sont ensuite regroupés en régions significatives 102 et 15 les pics présents dans lesdites régions sont ensuite analysés 103, les résultats étant présentés en sortie, par exemple sur un écran. Des exemples de mise en oeuvre de ces traitements sont détaillés plus loin dans la description à l'aide des figures 2 à 4. 20 La figure 2 donne un exemple de méthode de construction et détection des canaux significatifs. La méthode décrite ci-après requiert l'utilisation des distributions a posteriori obtenues par la méthode d'inférence Bayésienne semi paramétrique, c'est-à-dire par l'échantillonneur de Gibbs décrit 25 précédemment. Les canaux d'accumulation sont traités les uns après les autres, le canal d'indice j noté Bi étant le canal courant. Après échantillonnage de Gibbs, l'amplitude p des masses de Dirac, la position Z desdites masses de Dirac et l'amplitude q du fond sont utilisées par la méthode de construction et d'analyse de spectre 30 échantillonné selon l'invention. Après initialisation 200, un spectre de pics échantillonné en énergie par canal d'accumulation est construit 201 pour permettre la détection de régions significatives comprenant des pics. Le principe de détection des régions significatives repose tout 35 d'abord sur un test statistique effectué sur chaque intervalle d'énergie B; constituant un canal. Puis, on définit une région significative comme un ensemble de canaux contigus pour lesquels le test statistique est positif. Le nombre L de tirages Monte Carlo de l'échantillonneur de Gibbs définit la taille de la population des événements utilisés pour réaliser le test statistique. Une valeur de L limitée ne garantit pas de retrouver des masses de probabilité issues des tirages du processus de Dirichlet en nombre important sur chaque intervalle B. d'une région significative. Une condition garantissant que le nombre de tirages Monte Carlo s'avère suffisant est que ~» T_ où T. est le nombre de canaux de la plus grande région 1 o significative. En pratique, l'utilisateur ajuste le nombre L en fonction de la taille des régions significatives obtenues. Pour chaque tirage aléatoire de la procédure MCMC et pour tout j , on définit la quantité I -1j comme étant la masse de probabilité du canal 15 B ; du spectre normalisé de pics : niPk1(ZkEBJ) (1) k =1 expression dans laquelle : 20 j est l'indice associé à un intervalle B. donné ; N est le nombre de composantes (masses de Dirac) du processus de Dirichlet pour le tirage considéré; Zk est la position de la composante k du processus de Dirichlet produit par la procédure MCMC ; 25 la fonction 1(Cond) = 1 si la condition Cond est vraie et 1(Cond) = 0 sinon. En parallèle, on définit cDj la masse de probabilité du canal B . du spectre normalisé de fond pour chaque j par : 30 (2) Le spectre normalisé de pics peut ensuite être estimé 201, de manière à ce que celui-ci puisse être utilisé pour la détection de régions significatives comprenant des pics. L tirages issus de l'échantillonneur de Gibbs sont disponibles. Dans la suite de la description, la notation III) représente les valeurs générées de la variable aléatoire III à l'itération 1 de l'échantillonneur. De même, la notation 43;» représente les valeurs générées de la variable aléatoire 4^j à l'itération 1 de l'échantillonneur. Un estimateur du spectre de pics est choisi, par exemple, comme 1 o étant la moyenne conditionnelle de III pondérée par le nombre moyen de photons enregistrés affectés aux pics (w. n) où n est le nombre de photons enregistré et wE [0,1] est la probabilité du spectre de pics dans le spectre complet défini comme étant la somme pondérée du spectre normalisé de pics et du spectre normalisé de fond. A l'exemple d'application choisi est 15 celui de la spectrométrie Gamma, mais la méthode décrite peut être utilisée dans le cadre de la spectrométrie X ou de la spectrométrie de masse à temps de vol, par exemples. Pour tout canal d'accumulation d'indice j , des estimateurs Monte- Carlo sont construits en utilisant l'expression suivante : 20 L E(wnII1 IC)= nlw L ~-~ dans laquelle : C représente l'histogramme observé contenant dans chaque canal B. la 25 somme des événements (photons) dont l'énergie appartient à B1 . C peut être aussi désigné par les expressions « données » ou « observations » ; Pour ce qui est de l'estimation du fond, l'expression suivante est utilisée : E((1-w)n(D 1 1 C) = n 1(1- L dans laquelle : (3) 30 (1) (4) w l' représente le poids du spectre de pics pour le l ième tirage du processus MCMC. Les valeurs de E(w n II 1 1 C) et de E((1ù w) n cl) 1 1 C) peuvent, par exemple, être présentées graphiquement à l'aide d'un dispositif d'affichage. Ces estimateurs offrent à l'utilisateur une représentation globale de la densité de la distribution en énergie déconvoluée. Avantageusement, la méthode proposée permet, au-delà du seul calcul de la moyenne a posteriori, de déterminer un intervalle de crédibilité, par exemple à 95%, sur chaque intervalle Bi. Ainsi, l'incertitude associée à chaque canal énergétique est estimée, par exemple en définissant un intervalle de crédibilité à 95% [1_,îi+] déterminé en utilisant les expressions : = min(Pr(w n II 1 _ 0,025) - n ri (1'1' (LO,025L]) (5) F1+ = min(Pr(w n n 1 _ 0,975) - n{wW» fl }î (L0,975L]) (6) dans lesquelles : • pour tout x , la notation lx}" représente la collection des échantillons x triée par ordre croissant ; • pour tout x et pour tout j la notation lx}"( j) représente le j ième échantillon de la collection x triée par ordre croissant. • Ainsi, {w(»nr1T représente la collection des échantillons générés w(" f 1 triés par ordre croissant ; • L.] indique la partie entière ; • n représente le nombre total de photons enregistrés ; • r7 représente une variable caractérisant l'intensité du spectre de pic pour le canal d'énergie considéré. Une fois que le spectre de pics échantillonné en énergie est déterminé, un test de Kolmogorov-Smirnov 210 peut être appliqué. L'objectif de ce test 210 est de déterminer si sur chaque intervalle B . la distribution de la quantité E; correspondant à la probabilité du canal B; dans le spectre complet normalisé, et définie par l'expression : E, =wII;+(1ùw)cl) J (7) est significativement distincte de la distribution de F., F. étant une quantité correspondant au fond seul sur le j ième canal et étant défini par l'expression : I', _ (1-w)c (8) Pour ce faire, il est possible, par exemple, de mettre en oeuvre un test d'hypothèses non paramétrique permettant de déterminer si les deux échantillons Ei et F. associés au 1 ième canal d'accumulation suivent la même loi ou non. Le test de Kolmogorov-Smirnov à deux échantillons peut être utilisé dans le cadre de l'invention. Pour cela, la quantité DL(FJ,EJ) est déterminée en utilisant l'expression suivante : DL (T, , E,) = L ù sup 2 n CDF r.) } (q)ùCDF (1) 171 (9) dans laquelle : • {Ey } représente la collection des échantillons générés pour l sup B(CDF{r~ } (rJ) (10) expression dans laquelle B(.) est un pont Brownien, la distribution K = sup(B(4) étant définie par l'expression : x Pr(K x) =1ù 21(ù1) `- 2 Par conséquent, l'hypothèse que les deux échantillons E; et F. possèdent des distributions distinctes, c'est-à-dire que l'on est en présence d'un élément de pic significatif sur l'intervalle B; est vérifié si : DL(I'~,E.)> Ka (12) expression dans laquelle : 25 • Pr(K Le résultat du test de Kolmogorov-Smirnov détermine les étapes de traitements à appliquer ensuite. En effet : 20 - si DL(FJ,EJ )> Ka, il est considéré qu'un pic significatif a été détecté sur l'intervalle B; , et dans ce cas, des traitements 205 sont exécutés afin d'analyser ledit pic ; - si, au contraire, DL (FJ , E J) analyser 207. Si l'intervalle courant Bi ne correspond pas au dernier canal 3o d'accumulation à traiter, j est incrémenté 208 et un test de Kolmogorov- Smirnov est appliqué sur le nouvel intervalle courant. Si l'intervalle courant Bi correspond au dernier canal d'accumulation, le processus de traitement se termine. La figure 3 donne un exemple de mise en oeuvre du test de Kolmogorov-Smirnov dans le cadre d'une application de spectrométrie Bayésienne semi paramétrique en un canal Bi . Cet exemple comprend trois phases. La première phase est une phase d'initialisation 301, la seconde 311 correspond à une boucle de calcul et la troisième 312 à une phase dite de finalisation. La phase d'initialisation a pour but d'initialiser des variables intermédiaires utilisées par les deuxième 311 et troisième 312 phases. Ces variables sont notamment les variables lr, lE et dmaxet sont initialisées 301 telles que : lr = 0 lE =0 d,nax = 0 • IF correspond à l'indice de parcours de la collection CDF CL) (lr); [r. • lE correspond à l'indice de parcours de la collection CDF Os); {E, • dmax = max(lr ûix) correspond au décalage maximum entre les deux indices précédemment définis. La seconde phase 311, c'est-à-dire la boucle de calcul, est ensuite exécutée. Dans la suite de la description, pour tout intervalle Bi , la collection {Ey}T correspond à la collection d'échantillons {E(;) } obtenue par échantillonnage de Gibbs et ordonnée par ordre croissant. En suivant le même principe de notation, pour tout intervalle Bi , la collection {Fr }T correspond à la collection d'échantillons {Fr } obtenue par échantillonnage de Gibbs et ordonnée par ordre croissant. Un premier test 302 est ensuite appliqué de manière à vérifier l'inégalité suivante : {F? }î (lr) {E(1) }Î (lE ) (14) Où {I';') }T (lr) correspond au lr 1ème échantillon de la collection {Fr} triée par ordre croissant. Par analogie, Ey}T (ls) est le l~'ème échantillon de la collection {E») } triée par ordre croissant. Si l'inégalité (14) est vérifiée, la variable IF est incrémenté 303, c'est-à-dire que IF =lr +1. Si l'inégalité (14) n'est pas vérifiée, la variable ls est incrémenté 304, c'est-à-dire que ls = ls +1. 1 o Un second test 305 est ensuite appliqué de manière à vérifier l'inégalité suivante : dm 15 Si l'inégalité (15) est vérifiée, la variable dmax est mise à jour 306 telle que, dmax = lr ù ls . Puis un troisième test 307 est appliqué. Si l'inégalité (15) n'est pas vérifiée le troisième test 307 est appliqué directement. Le troisième test 307 a pour but de vérifier si l'inégalité suivante 20 est vérifiée : lr 1 d >Ka (17) Une variable de sortie vj booléenne est définie telle que celle-ci prend la valeur « 1 » si la présence d'un élément de pic significatif a été détectée sur l'intervalle B B. et la valeur « 0 » sinon. La valeur de vj est déterminée en fonction du résultat du test 308 vérifiant l'inégalité (17). Si l'inégalité (17) est vérifiée 310, v. =1 et la présence d'un pic est détectée sur l'intervalle B .. Si l'inégalité (17) n'est pas vérifiée 309, v.=0. 1 o La figure 4 donne un exemple de traitement de construction séquentielle des régions de pics significatifs. Une région de pics significatifs correspond à une région comprenant au moins un pic significatif, ladite région étant déterminée par détermination des canaux d'accumulation B B. contigus pour lesquels v; =1. 15 Cette opération correspond à une segmentation des canaux énergétiques. L'objectif de ce traitement est de fournir une liste de M régions R(m) (mE [1,...,M]), définies par leurs indices minimum jm et maximum j,, , lesdits indices correspondant à des indices de canaux d'accumulation. Les régions R(m) peuvent être déterminées en utilisant l'expression suivante : 20 3m R(m) = U B 1 Les indices minimum jm et maximum j, ainsi que le nombre de régions identifiées M sont déterminés par l'exécution des traitements décrits 25 ci-après. Une première étape 400 a pour but d'initialiser des variables intermédiaires. Ces variables sont notamment les variables j correspondant à l'indice de canal d'accumulation, m et 8 qui représentent respectivement l'indice de la région significative courante et la taille en nombre de canaux de 30 ladite région. Ces trois variables sont mises à zéro. Un premier test 401 est ensuite appliqué de manière à vérifier si v. =1. En d'autres termes, il est vérifié que le canal d'accumulation courant comporte un pic significatif ou une portion de pic significatif. Si v. =1, la variable 8 est incrémentée 402, c'est-à-dire 8 = 8+1. Sinon un second test 403 est appliqué, ledit test ayant pour objectif de vérifier si v._1 = 1, c'est-à-dire de vérifier si le canal d'accumulation précédent le canal courant contient un pic ou une portion de pic significatif. Si vj_1 =1, les variables jm , j,r , 8 et m sont mises à jours 405 de la manière suivante : 15 Un troisième test 406 est ensuite appliqué, ledit test vérifiant l'inégalité suivante : j L'estimateur d'intensité pour la région m est alors donné par la moyenne a posteriori : 2= ùD( (19) m L 1 m Un intervalle de crédibilité à 95% pour 2m est donné par les quantiles de la collection {42} en utilisant, par exemple, l'expression suivante : CI95% (2m) ~{2m~ }T (L0,025L]),{2m) } (LO,975L])] (20) dans laquelle la collection {2m`' }T correspond à la collection d'échantillons {2 } ordonnée par ordre croissant. La technique Bayésienne semi paramétrique permet de calculer directement tout autre moment ou fonctionnelle de la distribution de 2m . Il est également possible d'estimer 501 la localisation énergétique du centroïde de chaque région de pics significatifs. Pour cela, une estimation y)-k(» de la moyenne de la distribution l'énergie restreinte à Pm est déterminée. Pour tout m N = 11(Zk1) E Rm k =1 Pour tout m L Im L I1(Nm~) > o) l=1 Le dénominateur de l'expression (24) correspond aux nombre de 25 tirages MCMC pour lesquels il y a au moins une composante appartenant à la région Rm , les autres tirages ne pouvant pas être pris en compte dans l'estimation Monte Carlo. Rm , pk (1) peut être (21) 10 15 (22) l=1 (24) L'intervalle de crédibilité à 95% pour ,um est donné par les quantiles de la collection {,u,,,(»} et peut être déterminé, par exemple, en utilisant l'expression suivante : Cl95% (,um) = [{,u » }T (L0,025L]),{,uä,» }T (LO,975L])j (25) dans laquelle la collection {,uä,` }T correspond à la collection d'échantillons {,u,,,(» } ordonnée par ordre croissant. La technique Bayésienne semi paramétrique permet de calculer 1 o directement tout autre moment ou fonctionnelle de la distribution de Les régions de pics peuvent aussi être analysées en estimant l'étendue de chaque région de pics significatifs 502. Plusieurs quantités peuvent être fournies afin d'appréhender cette 15 grandeur. Par exemple, une estimation de l'écart-type am de la distribution d'énergie restreinte à la région Rm peut être utilisée. Pour tout m La technique Bayésienne semi paramétrique permet de calculer directement tout autre moment ou fonctionnelle de la distribution de am . L'écart-type am, dans le cas d'une loi gaussienne sur Rm peut être un indicateur de l'étendue de Rm notée ÔRm qui peut être déterminée en 1 o utilisant par exemple l'expression suivante : 56 (29) Étant donné que la restriction de la distribution d'énergie à tout Rm 15 peut être nettement différente d'une loi gaussienne, l'estimateur 6m peut ne pas être pertinent pour estimer l'étendue de la région. Une estimation de l'intervalle de plus forte probabilité à 99% de la restriction à Rm de la densité d'énergie, notée Am,99% , peut être estimée. Am,99% correspond à la largeur de la région de pics prise en sa 20 base. Pour cette quantité un estimateur basé sur sa moyenne conditionnelle et un intervalle de crédibilité 95% sont proposés. Pour tout m procédure MCMC est estimée la fonction cumulative de la distribution de l'énergie CDFI(» O restreinte à Rm en utilisant par exemple l'expression : 25 N CDFRim ( ) = 1)1(Zki~ aE [0,1], Q(Ri)(a)=fia ~CDF»( )=a. L'estimateur Ôm,99% de Am,99% , équivalent à l'intervalle de plus grande 5 probabilité à 99% de la distribution d'énergie restreinte à Rm , est alors donné par la moyenne a posteriori de O(m,99% , c'est-à-dire par l'expression : L IA(1) m,99% 1=1 1 o L'intervalle de crédibilité à 95% pour Am,99% est donné par les quantiles de la collection {A(,,,»,99% 1, c'est-à-dire par l'expression : CI95% (Am999% ) = [{gm),99% }î (LO,025L]),{O(m),99% }' (LO,975L])] (33) 15 dans laquelle la collection {O(m)99% }T correspond à la collection d'échantillons {0(ä ),99% } ordonnée par ordre croissant. La technique Bayésienne semi paramétrique permet avantageusement de calculer directement tout autre moment ou fonctionnelle de la distribution de Am,99% . I1(Nm~) > 0) l=1 (32) 20 REVENDICATIONS1- Procédé d'analyse de mesures spectrométriques, lesdites mesures C = {cl,...,cY} étant classées dans un histogramme, ledit histogramme étant composé de canaux d'accumulation, un canal j correspondant à un intervalle d'énergie Bi , ledit procédé étant caractérisé en ce qu'il comporte au moins trois étapes de traitement : - une première étape (100) déterminant par canal d'accumulation j des distributions p représentatives de l'amplitude et Z représentatives de la position des impulsions de Dirac composant le spectre normalisé de pics ainsi que des distributions q représentatives de l'arbre de Pôlya caractérisant le spectre normalisé de fond, lesdits éléments étant obtenus par échantillonnage de Gibbs ; - une seconde étape (101) de détection des canaux significatifs, ladite détection étant réalisée par application d'un test de Kolmogorov-Smirnov par canal d'accumulation de l'histogramme sur la base des résultats de la première étape (100) ; - une troisième étape (102) d'identification de régions significatives du spectre par regroupement des intervalles Bi des canaux significatifs identifiés lors de la seconde étape (101), lesdites régions comportant au moins un pic significatif. 2- Procédé selon la revendication 1 caractérisé en ce que le spectre de pics est déterminé (202) pour tout canal d'accumulation d'indice j , en utilisant des estimateurs Monte-Carlo définis par une expression telle que : L E(w n II 1 I C) = n w(1) II (i) L ~-~ dans laquelle : II;'e correspond à la masse de probabilité du canal B . pour le ième tirage du processus MCMC ;L correspond au nombre de tirages Monte Carlo de l'échantillonneur de Gibbs ; wE [0,1] est la probabilité du spectre de pics dans le spectre complet défini comme étant la somme pondérée du spectre normalisé de pics et du spectre normalisé de fond ; w l' représente le poids du spectre de pics pour le l ième tirage du processus MCMC ; n est le nombre de photons enregistré. 3- Procédé selon l'une quelconque de revendications précédentes caractérisé en ce que le spectre de fond est déterminé pour tout canal d'accumulation d'indice j , en utilisant l'expression suivante : E((1ùw)n(D 1 C) = n 1(1- L dans laquelle : 4 1) représente les valeurs générées de la variable aléatoire cl)i à l'itération 1 de l'échantillonneur, J représentant la masse de probabilité du canal B; du spectre normalisé de fond pour chaque j 20 4- Procédé selon l'une quelconque de revendications précédentes caractérisé en ce que l'incertitude associée à chaque canal énergétique est estimée en définissant un intervalle de crédibilité à 95% [1_,i+] déterminé en utilisant les expressions suivantes : 25 = min(Pr(w n II 0,025) - n{w(Il (fi) }T (L0,025L]) =min(Pr(wnIIi 5- Procédé selon l'une quelconque de revendications précédentes caractérisé en ce qu'une quantité DL(FJ,EJ) est déterminée pour application du test de Kolmogorov-Smirnov (210) en utilisant l'expression : DL(r,,E)= L ù sup 2 n CDF{r(~) } (1/) ûCDF{E(1) } (17 dans laquelle : E; correspond à la probabilité du canal Bi dans le spectre complet normalisé ; F. est une quantité correspondant au fond seul ; {E~~~ } représente la collection des échantillons générés EY) pour 1 L ; 20 {Fis" } représente la collection des échantillons générés Fi" pour 1 L ; CDF V)})la fonction cumulative empirique de la distribution de E telle que CDFEY)}(17)= 1 I(: ) 6- Procédé selon la revendication 5 caractérisé en ce que la présence d'un élément de pic significatif sur un intervalle B . est détectée si DL(Fi,Ej )> Ka, Ka étant défini tel que Pr(K Ka) =1-a, a correspondant à un niveau choisi de signifiance. 7- Procédé selon l'une quelconque des revendications précédentes caractérisé en ce qu'il comporte une étape d'estimation de l'intensité des pics sur une région significative Rm , ladite intensité étant déterminée en utilisant les expressions : L 2 m = -12(1) mavec 2( ) = n w(i) l'',, i) 1(z," E Rm L ld kd dans lesquelles : la fonction 1(Cond) =1 si la condition Cond est vraie et 1(Cond) = 0 sinon ; {R,,,} la collection de régions significatives du spectre ; 15 Zk est la position de la composante k du processus de Dirichlet produit par la procédure MCMC ; 8- Procédé selon la revendication 7 caractérisé en ce qu'il comporte une 20 étape de calcul d'un intervalle de crédibilité à 95% pour 2 en utilisant l'expression : C195% (2m) [{2) }T (L0,025L]),{2m) }T (LO,9 / 5L])] 25 dans laquelle la collection {2,`' }T correspond à la collection d'échantillons {2.m' } ordonnée par ordre croissant. 9- Procédé selon l'une quelconque des revendications précédentes caractérisé en ce qu'il comporte une étape d'estimation du centroïde 30 ,ûm de chaque région de pics significatifs en utilisant l'expression suivante :L E(l) _ l=1 P. L Il(Nm(l) 1=1 dans laquelle : ,u,,,(» est une quantité estimée sur chaque région significative Rm et pour chaque tirage de la procédure MCMC en utilisant l'expression : N /e - 1pk O zk"I(zk" E Rm ; k =l le dénominateur correspond aux nombre de tirages MCMC pour lesquels il y a au moins une composante appartenant à la région Rm , les autres tirages ne pouvant pas être pris en compte dans l'estimation 1 o Monte Carlo. 10- Procédé selon la revendication 9 caractérisé en ce qu'il comporte une étape de calcul d'un intervalle de crédibilité à 95% pour ,um en utilisant l'expression : 15 CI95%(Jlm)= L{ }T (L0,O25LJ),{,um~> }T (Lo,975L])j dans laquelle la collection {du }T correspond à la collection d'échantillons {,u,n(» } ordonnée par ordre croissant. 11- Procédé selon l'une quelconque des revendications précédentes caractérisé en ce qu'il comporte une étape d'estimation dm de l'écart- type am de la distribution d'énergie restreinte à la région Rm en utilisant l'expression : 25 L lam(1) l=1 6m _ L ~1~Nmi) o) 1=1 dans laquelle : 33 20(Z »21(e E Rm ) û N k \ k=1 1(Zk" E Rm ) 2 N 12- Procédé selon la revendication 11 caractérisé en ce qu'il comporte une étape de calcul d'un intervalle de crédibilité à 95% pour am en utilisant l'expression : CI95% (cm) = [{6ä,l) }l (L0,025L]),{6mw }l (L0,975L])] dans laquelle la collection {6ä~» }T correspond à la collection 1 o d'échantillons {6m(1' } ordonnée par ordre croissant 13- Procédé selon l'une quelconque des revendications précédentes caractérisé en ce qu'il comporte une étape d'estimation de la largeur de la région déconvoluée en sa base, définie comme l'intervalle de 15 plus forte probabilité à 99% de la restriction à Rm de la densité d'énergie (notée Am,99%) , en utilisant l'expression L IA(1) m,99% Ô = l=1 m,99% L 11(N,, l=1 20 dans laquelle : A(1) 9% = [Q,(:2 (0,005),Q(RI) (0,995)] où Q(RI, (a) = CDFI, ( ) = a et CDF " ( ) = D3 l) 1(Zk1) 14-Procédé selon la revendication 13 caractérisé en ce qu'il comporte une étape de calcul d'un intervalle de crédibilité à 95% pour Am,99% en utilisant l'expression : CI95% (Am,99%) = [{A ,99% 025L]),{gm",99% }«L0,9751-1dans laquelle la collection {A 99% }T correspond à la collection d'échantillons {0(äi'99%} ordonnée par ordre croissant. 15- Dispositif d'analyse de mesures spectrométriques, lesdites mesures C = {cl,...,cY} étant classées dans un histogramme, ledit histogramme étant composé de canaux d'accumulation, un canal j correspondant à un intervalle d'énergie B1 , ledit dispositif étant caractérisé en ce qu'il comporte des moyens pour mettre en oeuvre le procédé selon l'une 1 o quelconque de revendications 1 à 14.