FR3034221A1 - Procede et dispositif pour detecter des radioelements - Google Patents

Abstract

L'invention concerne un procédé pour déterminer la nature des radioéléments présents dans un objet et leur activité caractérisé en ce qu'il comporte au moins les étapes suivantes : • une première phase (31) de simulation numérique de réponses spectrométriques pour un ensemble d'énergie incidente E et un ensemble d'énergie de sortie mesurée E', afin d'obtenir un ensemble de données simulées, • une deuxième phase de régression non paramétrique (32) sur les données simulées, estimation non paramétrique de la quantité représentant la probabilité jointe des triplets (E,E',y) à partir de points simulés (Ei,Ei',yij) afin d'en déduire un méta-modèle S(E, E') pour tout couple d'énergie (E, E') sur une fonction continue, • à partir du méta-modèle S(E, E'), (34) la détermination de la nature et de l'activité des radioéléments présents dans l'objet.

Description

1 PROCEDE ET DISPOSITIF POUR DETECTER DES RADIOELEMENTS L'objet de l'invention concerne un procédé et un dispositif pour détecter des radioéléments contenus dans un objet, déterminer leur nature et leur activité. Elle permet notamment de tenir compte d'effets parasites provenant de l'environnement dans lequel se trouve l'objet à analyser. La spectrométrie gamma est une technique permettant de détecter les rayonnements gamma émis spontanément par des radioéléments. Cette technique permet d'identifier les radioéléments émetteurs gamma et de mesurer leur activité. Pour atteindre des objectifs qualitatif et quantitatif, le traitement d'un spectre gamma se base sur l'analyse des pics présents dans le spectre. Ces derniers correspondent au rayonnement gamma émis par un radioélément spécifique, par exemple, le 137Cs émettant un gamma à 661,5 keV ou le 60Co émettant un gamma à 1173,2 keV et un gamma à 1332,5 keV. La présence d'un pic dans un spectre gamma signifie que l'énergie du rayonnement gamma incident a été entièrement absorbée dans le détecteur, suite à une ou plusieurs interactions, effet photoélectrique, effet Compton ou production de paires, ayant eu lieu au sein de ce dernier.

Les détecteurs utilisés pour réaliser des mesures de spectrométrie gamma peuvent être séparés en deux catégories : les scintillateurs et les semi-conducteurs. La qualité d'un spectre gamma peut être évaluée grâce à deux paramètres propres au détecteur utilisé: la résolution en énergie, c'est-à-dire la capacité à séparer deux pics proches en énergie définie comme étant la largeur à mi-hauteur d'un pic présent dans le spectre à une énergie donnée, et l'efficacité de détection, plus le matériau sera dense, plus la probabilité d'absorber l'intégralité du rayonnement gamma incident est importante. Ces deux paramètres varient de manière significative en fonction du détecteur utilisé.

Les scintillateurs inorganiques (Nal, BGO) bénéficient d'une résolution dégradée par rapport aux détecteurs semi-conducteurs (CZT ou 3034221 2 HPGe). Néanmoins, ils présentent généralement une plus grande efficacité de détection. En effet, en raison de leur faible coût, il est possible de les fabriquer dans de grandes dimensions. La figure 1 montre, dans un diagramme Energie de canal, coup par canal, un exemple de spectre S(I) 5 obtenu à l'aide d'un détecteur HPGe et d'une source de 137Cs émettant un rayonnement gamma à 661,5 keV. La figure 2 compare des spectres S(HPGe), S(CZT), S(Nal) obtenus à l'aide d'un échantillon d'uranium et de différents types de détecteurs, présentant des performances différentes en termes de résolution en énergie.

1 0 Il est aussi connu d'utiliser des scintillateurs plastiques pour détecter des rayonnements gamma. Le coût de fabrication de tels scintillateurs étant assez faible, ces derniers peuvent être fabriqués dans de grandes dimensions pour différents types d'application. Du fait de leurs caractéristiques intrinsèques et en dépit de leurs dimensions très importantes 15 par rapport à d'autres détecteurs, la probabilité qu'un photon incident dépose toute son énergie au sein d'un détecteur par effet photoélectrique est extrêmement faible. La figure 3 montre un spectre S(II) obtenu à l'aide d'un scintillateur plastique de type EJ200, en présence de sources 60Co et 137Cs. Les formes caractéristiques de ces spectres sont dues aux fronts Compton 20 faisant suite aux interactions au sein du détecteur. On peut ainsi noter l'absence des pics photoélectriques caractéristiques de la présence de 137Cs et de 60Co. Pour cette raison, ces détecteurs sont généralement utilisés pour réaliser des mesures de comptage total, détection de tous les rayonnements gamma interagissant au sein du détecteur sans information sur leur énergie.

25 Une approche traditionnelle consistant à mesurer l'aire nette des différents pics d'intérêt n'est donc pas envisageable et une solution alternative doit être employée si l'on souhaite procéder à l'identification des radioéléments présents dans le spectre. Une solution pour traiter un spectre gamma, identifier les 30 radioéléments émetteurs et mesurer l'activité de ces derniers, consiste à analyser le spectre dans son ensemble et à ne pas se restreindre 3034221 3 uniquement à l'analyse des pics photoélectriques. L'approche suivie consiste à exprimer l'analyse du spectre sous la forme matricielle suivante : S=H.A, avec S : une matrice signal, matrice de dimensions (nbre_canaux, 1), où la 5 variable nbre_canaux est égale au nombre de canaux du spectre gamma, A : une matrice activité, matrice de dimensions (nbre_energies_incidentes, 1). La variable nbre_energies_incidentes correspond au nombre d'énergies incidentes défini par l'utilisateur et pris en compte dans la reconstruction. Ces énergies incidentes correspondent au nombre d'éléments de volume 10 élémentaires d'un problème classique de tomographie d'émission, H : un projecteur du problème, matrice de dimensions (nbre_canaux, nbre_energies_incidentes). Cette matrice comprend l'ensemble des efficacités de détection intervenant dans le problème. A titre d'exemple, un élément hij de cette matrice correspond à la probabilité qu'un photon 15 d'énergie incidente j soit détecté dans le canal i du spectre gamma. La matrice signal S correspond au résultat de la mesure, au spectre gamma à analyser, la matrice projecteur est généralement calculée par simulation, la matrice A correspond au résultat de la reconstruction. La figure 4A et la figure 4B illustrent un exemple de notion de 20 projecteur H, pour une énergie incidente de 660keV, figure 4A. Le spectre S(IV) de la figure 4B correspond à l'énergie incidente correspondant à une colonne entière du projecteur H. Le projecteur contient donc la réponse spectrale du détecteur pour chaque énergie incidente prise en compte dans le problème et définie par l'utilisateur. C'est sur cette grille d'énergies 25 incidentes que sera effectuée l'étape de reconstruction. Différentes méthodes permettent de prendre en compte l'ensemble du spectre gamma afin de procéder à une analyse qualitative et quantitative, permettant l'identification des radioéléments présents dans le spectre et l'activité du radioélément présent dans le spectre des énergies 30 incidentes. Cette méthodologie d'analyse est illustrée sur les figures 5 et 6. La figure 5 représente le spectre gamma à analyser, correspondant à une 3034221 4 simulation de la réponse d'un scintillateur plastique de type EJ200, en présence de trois sources gamma (241Am, 137C s, Rn S --Co). Le spectre reconstruit Sr à partir des résultats de l'analyse et des énergies incidentes est représenté en traits pointillés sur la figure. La figure 6 illustre le résultat 5 d'une reconstruction effectuée à l'aide d'un algorithme de type ML-EM, les énergies incidentes reconstruites à partir du spectre de la figure 5. Les méthodes d'analyse traditionnelles malgré tout leur intérêt et avantages présentent quelques limitations. Une des limitations intrinsèques est liée au caractère discret de la grille choisie pour effectuer la 10 reconstruction et obtenir les énergies incidentes. Dans toutes les approches classiques, comme celle de type ML-EM, cette grille est discrète, suivant un pas défini par l'utilisateur. Ce pas est l'une des limitations lors de l'étape de reconstruction puisqu'il ne sera possible de reconstruire que sur la grille définie initialement. En fonction de la complexité du problème traité, le pas 15 entre chaque énergie incidente peut être adapté. Un pas fin permettra une plus grande finesse lors de la reconstruction, mais aura un impact direct sur le calcul du projecteur. Ce dernier étant généralement calculé par simulation, un pas plus fin se traduira par une augmentation du temps de calcul, ce dernier pouvant devenir prohibitif pour certaines applications.

20 Une alternative consisterait à chercher les solutions sous forme d'un spectre discret dont le support est un espace continu. Une telle approche pour être intégrable dans le cadre d'un algorithme ML-EM demanderait à calculer le projecteur à la volée par simulation Monte-Carlo pour chaque énergie incidente au cours des itérations de l'algorithme.

25 Comme les valeurs d'énergies incidentes évoluent au gré des itérations, des calculs de réponses sont requis à chaque itération sans que les valeurs précédemment calculées puissent être réutilisées. L'ordre de grandeur de calcul d'une telle approche s'avère prohibitif pour les calculateurs actuels. Dans la reconstruction de type ML-EM, aucun a priori n'est introduit lors de 30 l'étape de reconstruction. Pour la représentativité du projecteur, une des principales limitations de ce type d'approche est liée à la représentativité du 3034221 5 projecteur utilisé lors de la reconstruction. Ce dernier est généralement obtenu suite à une étape de simulation, à l'aide d'un modèle décrivant le détecteur utilisé et son environnement. Néanmoins, il est généralement difficile d'aboutir à une modélisation fine, connaissance de la géométrie du 5 détecteur, limitations du code de calcul qui ne modélisent qu'imparfaitement les interactions au sein du détecteur, absence de prise en compte de certains phénomènes physiques, comme la collection de photons visibles par exemple dans le cas des scintillateurs plastiques. Si la modélisation du détecteur et/ou de l'environnement est imparfaite, un biais plus ou moins 1 0 important sera présent lors de la reconstruction. Une autre approche, plutôt que considérer une grille d'énergies incidentes sur laquelle on effectuerait la reconstruction est d'utiliser une base de données. Dans ce cas, on ne simule plus une grille d'énergies incidentes, mais directement la signature d'un élément donné, par exemple pour le 60Co 15 on simulera directement la réponse du détecteur pour deux énergies incidentes, émises respectivement à 1173,2 keV et 1332,5 keV. Avec cette approche la reconstruction est effectuée directement sur une grille de radioéléments. Elle présente aussi l'avantage de converger plus rapidement que l'approche avec grille, car le nombre de composantes à prendre en 20 compte est plus faible. L'un des inconvénients de cette méthode est sa limitation à ne prendre en compte que les radioéléments initialement présents dans la base de données. Si l'analyse d'un spectre gamma fait intervenir un radioélément absent de la base de données, la reconstruction sera erronée.

25 Dans la suite de la description, on désigne sous l'expression « spectre d'entrée » un ensemble d'énergies incidentes définies par un utilisateur et sous l'expression « spectre de sortie » un ensemble d'énergies mesurées. Le mot méta-donné et le mot projecteur désignent un même objet. L'invention concerne un procédé pour déterminer la nature des 30 radioéléments présents dans un objet et leur activité caractérisé en ce qu'il comporte au moins les étapes suivantes : 3034221 6 - une première phase de simulation numérique de réponses spectrométriques pour un ensemble d'énergie incidente E et un ensemble d'énergie de sortie mesurée E', afin d'obtenir un ensemble de données simulées, 5 - une deuxième phase de régression non paramétrique sur les données simulées, estimation non paramétrique de la quantité représentant la probabilité jointe des triplets (E,E',y) à partir de points simulés (Ei,Ei',yij) afin d'en déduire un méta-modèle S(E, E') pour tout couple d'énergie (E, E') sur une fonction continue, 10 - à partir du méta-modèle S(E, E'), la détermination de la nature et de l'activité des radioéléments présents dans l'objet. Selon une variante de réalisation, le procédé comporte les étapes suivantes : - on utilise n un nombre de points dans la grille d'entrée, énergies E, et 15 n' un nombre de points dans la grille de sortie, énergies E', - pour i = 1, ...,n et j = 1, ...,n', on dispose des données caractéristiques des intensités spectrales calculées À pour une énergie d'entrée EL et une énergie de sortie E'1, - on utilise une méthode d'estimation non paramétrique de la quantité 20 f(E,E',y) représentant la densité de probabilité jointe des triplets (E,E',y) à partir des points simulés (Ei,E'pyii), et on déduit un modèle S(E,E') pour tout (E, E') E rz2 où I': est un espace continu : S(E, = IE (yIE,E') fll, y . f(y1E,E)dy = fR y f(E,E',y) dy fR f (E, y) dy où le symbole IE(y) représente l'espérance mathématique de la variable aléatoire y, y est déduite des intensités spectrales calculées x,;.

25 On peut déterminer les valeurs À en utilisant un logiciel Monte- Carlo, les données À étant considérées comme des réalisations issues 3034221 7 d'une distribution de Poisson dont l'intensité ILS est estimée au moyen d'une procédure de régression non paramétrique et on introduit yii = + E) OÙ 0 E 1, puis on approxime la distribution de probabilité de yii pour des valeurs 5 suffisamment grandes de IL> > 10 par une loi gaussienne de moyenne log ILS et de variance pour la loi jointe f y), on choisit un mélange par processus de Dirichlet (DPM) comme distribution a priori et on exprime la distribution aléatoire en une somme sur l'infini de foi, composantes de f , y), f(E, y) =wk fok(E,E, y) k=1 10 paramétrée par 61k le paramètre associé à la kième composante de G mesure aléatoire définie par = 8ok(') k=1 avec w1 = V1, Wk = - V1) tel que Vk-Beta(1, a) et (Su(-) représente la fonction de Dirac localisée en u et Beta(a, b), pour 0 < x < 1 avec = rr((aa)-Er(bb))xa 1(1_ x)b-1. fBetea,b)(X) 15 Les composantes de la loi jointe foi, sont, par exemple, exprimées à partir des valeurs suivantes : - 0k = (ii>k,ik,ipk,flk) avec rik = = lek, k), - Xk(E) le vecteur de régresseurs centré avec É = (E,E'), et /3k le vecteur de coefficients de régression, 20 - la matrice Ek, dépendant du paramètre tpk E {0,1}, comme suit : Ek = k (01 (1 -h si R,h-1 (Tk T,k) .Rek avec Rip = °) Si tpk = 0 et Rek = k. 1 2 1 lPk = 1, permettant de choisir entre une composante alignée sur les axes E et E' (Ipk = 0) et une composante oblique orientée par la droite E = (IP k = 1). 3034221 8 chaque composante fek s'exprime alors : fOk(E, E, Y) = g I Tik, -7'r (YIK - g k(É), fiTc.4(É)) où, .7«. 62) représente la loi gaussienne de kt et de variance a2, , E) la loi gaussienne bivariée de moyenne fi E rz2 et de matrice de covariance E, où 5 la loi a priori du paramètre uk est une gaussienne, la variance -ck est distribuée suivant une loi gamma-inverse, tpk suit une loi a priori de type Bernoulli et l'a priori pour les coefficients de régression coefficients /3k est une loi normale (gaussienne) multivariée de dimension IX(É)I, et en appliquant une règle de Bayes sur l'expression de ffE,E',y) on obtient 10 l'expression f (yIE , E') Wk j'r2(EI> J^r (y lig - gk e-efek(É)) k=i L=1 W/ -7\r2n 1174,10 et le modèle probabiliste S(E, E)= IE (yIE , E) ov k N2 (É 111:' lk) le' - 'g (É>) JVZ(Lifi 1, El) k k ^ à partir de ce modèle probabiliste et des données observées par simulation , on estime la loi a posteriori f(E,E',yIE1,E'1,Y11. --- .En. , n' ynn,) et l'espérance conditionnelle .§.(E, E') = IE (yIE , , E1, E v --- En, E' n' ynn, ) 15 afin de déterminer les éléments présents dans l'objet et leur activité. Pour le calcul de l'a posteriori, on utilisera, par exemple, un schéma d'approximation Monte-Carlo par chaîne de Markov (MCMC), pour toute itération (t) de la procédure MCMC, on génère une réponse spectrale débruitée S(E,E')(0, pour T générations, la distribution a posteriori de la 20 réponse spectrale est approximée par l'ensemble des tirages S(E,E')(t) pour t = 1, T, et la réponse estimée s'exprime : T -,>-- T 1S(E , Er) t=1 3034221 9 Le schéma d'approximation peut comporter une étape d'échantillonnage par tranche utilisant un nombre aléatoire fini K de composantes pour chaque itération et en ce qu'il comporte les étapes suivantes : 5 on introduit des variables latentes de classification Kip définies pour i = 1, ...,n et j = 1, ...,n', tel que Ki = k si (Ei,ri,yii) est distribué suivant la kième composante du mélange f(E,E1,y), on définit un modèle pour les paramètres du mélange, pour tout i < n, j < Kiji W1, W2, --- ^ Wk 8k (') k=1 EL, EjlKij, 01, 02, --- - E il/11(w 1Ki1) 1 0 avec it> ri] ° (11Kii' 'Ki]) 1Kij ReKiJ. T ReKi; Kii,Ei, 01, 02, - N (yii I PLI - gKii(Ei, e -RKii équivalent à la vraisemblance f(K11, , E, E'n,, y , à partir de ces distributions de probabilités et par application de la règle de Bayes, on calcule la densité de probabilité conditionnelle f(w1, W2, --- 611, 02, --- I K11, --- Knn', E1, E 1, - - - , En, en', Y11, - - - , Y nn') en utilisant un échantillonneur de Gibbs qui à chaque itération (t) génère 15 successivement les échantillons suivants, pour tout k, W1, W2, --- I K11, --- Knn' kI Kll, --- Knn', E - -- , En, E T xi 11-k» Pk, K11, - - - Knn', E1, E ff, -- - , En, E --- Knn', E1, E -- - , En, en' µk'ki k, k, k, K11, - - - , Knn', --- En, en' le xi k, IN, K11, --- Knn', E1, E i, - - - , En, E in', Y11, - - - , Ynn' Ynn'i w1, w2, - - - , 01, 02, - - - ) 3034221 10 et au niveau des nombres de points i dans la grille d'entrée d'énergie et de points j dans la grille de sortie, pour tout i < n, j < Yi), [t2, --- 2, --- T1, T2, --- T'1, T'2, --- ,1(31,1(32, --- 11)1, 1P2, --- VV1, VV2, --- à partir de T itérations de l'échantillonneur de Gibbs, on obtient une estimation du méta-modèle pour la spectrométrie (E, E') -1S (E,E')(t) t=i 5 Il est possible d'introduire une variable auxiliaire uij afin de ne générer qu'un nombre aléatoire fini K de composantes à l'itération (t) tout en évitant une troncature arbitraire du modèle. Le procédé peut comporter une étape de calcul de l'écart-type a posteriori et des intervalles crédibles à partir de l'ensemble {S(E,E')(t)} 1 65 1 (E,E) 1 (E, El) -S (E , e)(0)2 - t=1 10 Selon une variante de réalisation, on introduit un a priori Gamma(cpb,4,) sur les paramètres d'échelle br et b' dans la distribution des amplitudes de composantes conduisant à une distribution a posteriori Gamma : / Kn 1 br- Gamma (Pb Kn, k=1 avec Gamma(a, b), pour x > 0, fGammea,b)(X) - r(a) x e ba a-1 -bx 15 L'a priori gaussien sur les coefficients de régression /3k peut être remplacé par un a priori basé sur le tirage aléatoire de P points E' P) à partir de .7V'2(fik, Ek) et on prend pour (Yi, ,5"p) la plus proche valeur de yii correspondant à chaque point échantillonné où P est 3034221 11 plus grand que la taille du vecteur pk, on génère alors J^r(m,g, Fp) comme loi a priori (avec Ép = p)) r I= P \1 /3)el<gpir'gkgp .e-S1P \p=1 /P Mi& = rig - gk(Ép) e-519 \p=1 Selon une variante de réalisation, le procédé génère un méta-modèle étendu noté S(E,É,) où est un paramètre identifié par un indice 5 entier et caractéristique d'un effet de matrice. Le procédé peut estimer l'activité des radioéléments en exécutant les étapes suivantes : soit Nx le nombre d'émetteurs retenus pour l'élément x considéré et n-x,/ pour / = 1, ...,Nx, les probabilités d'émissions associées ainsi que vx,/ pour 10 1 = 1, ...,Nx, les énergies correspondantes, on définit la réponse du radionucléide tpx(E, , pour une énergie observée E' et un effet de matrice par Nx tp Trx,i S(v xj, E, 1=1 ir on définit = fo tpx.e') de et rcx*/ = pour tout / = 1, ...,Nx, x,e on définit une réponse de radionucléide normalisée par NX tpx* ,(e) =71-x S(v xj, E', 1=1 15 on vérifie que fo- tpx* (El) dr = 1 on exprime la densité de probabilité du ième photon observé dans le canal d'énergie E; pour i = 1, ..., n f(Ei)=wkiPXk,fic(P'E;) k=1 3034221 12 où p est un paramètre positif de conversion de l'indice de canal en énergie (keV/canal). On introduit les variables Ki d'allocation du -ème photon observé à une composante k du mélange, 5 On alterne des tirages aléatoires suivant des lois conditionnelles suivantes, pour tout i < n, pour tout k w 2, -1 Kn Xkl Kn,E;, , En, ---, P kI K1, Kn, E;, En' , Xi, X2, P ainsi que pl K1, , , En,X1,X2, ---, --- en utilisant un nombre aléatoire fini de composantes à l'itération (t), 10 à partir de T itérations de l'échantillonneur de Gibbs on obtient une estimation des activités radioéléments intervenant dans le mélange pour tout k, T 1 A- A(t) Xk T k Xk k t=1 avec p un a priori gaussien centré sur µP et d'écart-type ap. P J'r (P I Pp, ).

15 On calcule par exemple l'écart-type a posteriori des activités en exécutant les étapes suivantes : 1 T G 1- w(t)A(t), )2 AXk - 1 L_I x k k k/ t=1 et/ou le spectre d'entrée estimé, déconvolué de la réponse du système N (t) T Xk 1 (t) S .(E) k ITxk (t) 1 (t) 15vt) (E) X ,1 t=1 1=1 c )2 3034221 13 D'autres caractéristiques et avantages de la présente invention apparaîtront mieux à la lecture de la description qui suit donnée à titre illustratif et nullement limitatif annexée des figures qui représentent : - La figure 1, un exemple de spectre gamma obtenu avec un détecteur 5 HPGe, - La figure 2, un comparatif de spectres d'uranium obtenus à l'aide de différents types de détecteurs, - La figure 3, un exemple de spectre obtenu à l'aide d'un scintillateur plastique, 10 - La figure 4a et la figure 4B, un exemple d'énergie incidente à 660 keV et le spectre mesuré par un scintillateur plastique correspondant à cette énergie incidente, - La figure 5, un exemple de spectre gamma à analyser et un spectre reconstruit selon l'art antérieur, 15 - La figure 6, les énergies reconstruites à partir du spectre présenté sur la figure 5, - La figure 7, un exemple de système pour la mise en oeuvre du procédé selon l'invention, - La figure 8, un synoptique des étapes mises en oeuvre par le procédé 20 selon l'invention, - Les figures 9 et 10, un exemple de spectre des énergies incidentes obtenu par le procédé selon l'invention, - Les figures 11 et 12, un deuxième exemple de spectre des énergies incidentes obtenu par la mise en oeuvre de l'invention.

25 La figure 7 est un exemple de système permettant la mise en oeuvre du procédé selon l'invention. Le système comporte un objet 10 qui est susceptible d'émettre spontanément des rayonnements gamma du fait de la présence d'éléments radioactifs. Les rayonnements gamma sont reçus par un dispositif 20 comprenant un module de détection des rayonnements, 21, 30 relié à une électronique d'acquisition permettant d'obtenir un spectre, un processeur 23 ou unité informatique adapté à exécuter les étapes du 3034221 14 procédé selon l'invention, à traiter les données du spectre obtenu, un module d'affichage 24 des résultats, un moyen de stockage 25 pour la sauvegarde des résultats, par exemple. Le moyen de stockage peut aussi mémoriser les énergies incidentes définies par un utilisateur dans le contexte d'une 5 application donnée. Le procédé selon l'invention repose notamment sur une utilisation de l'ensemble de l'information présente dans le spectre à analyser afin de réaliser une succession d'étapes itératives de déconvolution du spectre incident. Le procédé utilise un projecteur ou méta-modèle continu. On va 10 rechercher la réponse du système de spectrométrie par une technique de simulation en calculant chaque fois qu'il est nécessaire la réponse complète en énergie du système pour une configuration donnée du spectre d'entrée. On va rechercher la réponse du système sous forme d'un spectre discret dont le support (énergies incidentes) est un espace continu. Ce modèle 15 continu va permettre de calculer la réponse du système de spectrométrie pour toute énergie incidente E et toute énergie de sortie observée E', méta-modèle noté S(E, E'). Une étape de reconstruction bayésienne utilise le méta-modèle pour déterminer les radioéléments contenus dans l'objet à analyser.

20 Le procédé va comporter une première phase au cours de laquelle on va effectuer une simulation numérique de réponses spectrométriques S, 32, pour une grille d'énergies incidentes Ei, 31, suivie d'une deuxième phase où l'on applique une régression non paramétrique (prédiction statistique), 33, sur les données simulées. Lors de la deuxième phase, on injectera, par 25 exemple de la connaissance physique, 34, sous forme d'a priori afin de guider la régression dans des régions de l'énergie incidente où on ne dispose pas de données simulées. En substance, ceci permet d'indiquer que sur des domaines restreints d'énergie, la signature de la réponse spectrale présente soit des caractéristiques proportionnelles à l'énergie d'entrée, soit 30 constantes quelle que soit l'énergie d'entrée Ei. A partir du méta-modèle 35, et de la connaissance a priori, il sera possible de remonter, par application 3034221 15 d'un algorithme de déconvolution, 36, à la composition du spectre 37. Le modèle de régression choisi étant un modèle non paramétrique, le procédé ne présume pas d'un quelconque modèle linéaire ou polynomial et pourra ainsi s'adapter à n'importe quelle non linéarité dans la réponse.

5 Le procédé va chercher à estimer une fonction continue de 12 en E et E' à partir d'un nombre potentiellement infini de paramètres, ou composantes du modèle. Pour mettre en oeuvre le procédé on utilise, par exemple, n le nombre de points dans la grille d'entrée (énergies E) et n' le nombre de 10 points dans la grille de sortie (énergies E'). Pour i = 1, ..., n et j = 1, ..., n', on dispose des données caractéristiques des intensités spectrales calculées À pour une énergie d'entrée EL (données définies par l'utilisateur) et une énergie de sortie E'1 (données mesurées). Pour illustrer le procédé selon l'invention, on suppose que les valeurs À ont été obtenues par 15 l'intermédiaire d'un logiciel de simulation Monte-Carlo. Par conséquent, les données des intensités spectrales À sont considérées comme des réalisations issues d'une distribution de Poisson dont l'intensité ILS sera estimée au moyen d'une procédure de régression non paramétrique. Comme Àii est alors une quantité positive ou nulle, on introduit yii + E) où 20 0 < E « 1. Dans ces conditions, la distribution de probabilité de yii peut être approximée pour des valeurs suffisamment grandes de ILS (typiquement > 10) par une loi gaussienne de moyenne log ILS et de variance L . 1i; La distribution de poisson sera exprimée de la manière suivante : Distribution PoissonM, pour n E N, Ân. fPoisson(Â) (ri) = e n 25 Le procédé selon l'invention repose sur l'estimation non paramétrique de la quantité f (E, E, y), représentant la densité de probabilité jointe des triplets (E,E1,y) à partir des points simulés (Ei,E'pyii), afin d'en déduire le modèle S(E, El) pour tout (E, E') E où rz est un espace continu : 3034221 16 S(E,É) = IE (yIE,É) fll, y. f(y1E,É)dy = fR y . ffE,É,y) dy fR f(E,E',y) dy où le symbole IE(y) représente l'espérance mathématique de la variable aléatoire y, ce symbole représente une espérance statistique utilisé plus loin dans la description pour caractériser l'intensité d'une énergie dans le spectre. Pour la loi jointe ffE,E',y), on choisit, par exemple, un mélange par 5 processus de Dirichlet (DPM) comme distribution a priori. Cette modélisation s'exprime à partir d'une distribution aléatoire G telle que G - DP(a, G0) où le symbole "-" signifie "est distribué suivant" et DP(a, G0) représente la distribution d'un processus de Dirichlet (Hjort N L, Holmes C, Müller P, Walker S G, Ghosal S, Lijoi A, Prünster I, The Y W, Jordan M I, Griffin J, 10 Dunson D B and Quintana F 2010 Bayesian Nonparametrics Cambridge Series in Statistical and Probabilistic Mathematics) avec une mesure moyenne Go et un paramètre a dit de concentration. Cette distribution de probabilité sur des distributions de probabilité aléatoires joue un rôle central en modélisation bayésienne non paramétrique. Une représentation courante 15 du DP est issue de Sethuraman (1994) qui exprime la mesure aléatoire G(-) comme : = Wk 8ok(-) k=1 avec w1 = V1, Wk = - V1) tel que Vk-Beta(1, a) et (Su(-) représente la fonction de Dirac localisée en u. Ici, Ok-GO représentent les paramètres associés à la kième composante de G. Le caractère non paramétrique est 20 issu du nombre de composante potentiellement infini intervenant dans la somme caractérisant l'a priori sur G. La distribution Dirichlet(a a2, --- CeK) est définie avec x = xl()' E I , xi > 0 pour 1 < i < K, xi < 1, xK = 1 - xi et ai> 0 pour 1 < < K, 3034221 17 nai+a2+----FaK) xal-1xa2-1 K-1 fDirichlet(ai,a2 ..... aK)(X -\ ) = .1 2 nai)na2)...naK) La distribution Beta(a, b), est définie pour 0 < x < 1 : r(a + b) a-1(1 - X)b-1. f(E, y) = wk y) k=1 où foi, indique la kième composante de f , , y), paramétrée par Bk.

5 Les composantes foi, dans le cas de l'application du DPGLM à la caractérisation du méta-modèle de réponse énergétique, sont écrites en introduisant les notations additionnelles suivantes : - 0k = (ilk,Tk,ipk,f3k) avec fik = = (rk, Tk), - Xk(E) le vecteur de régresseurs centré avec É = (E,E'), et /3k le vecteur 10 de coefficients de régression. On remarquera que le modèle de régression pour l'issue y sachant les covariables (E,E') peut être choisi linéaire (i.e. ge) = (1,E - - µ k)T où vT représente la transposée du vecteur y) ou polynomial dans l'approche DPGLM (ex. gk(É) = (1,E - - ktk)(E ' µ'k), (E 11-02 111JY pour une régression 15 quadratique), - De plus, on définit la matrice Ek, dépendant du paramètre tpk E {0,1}, comme suit : Ek = R k. /, T ,). Ripk avec Rej, = (01 °i) si tpk = 0 et k n ek -\77 1 1) !l = 2 1 (\ 1 1 si tpk = 1. La variable discrète binaire tpk permet ainsi de choisir entre une composante alignée sur les axes E et E' (Ipk = 0) et une 20 composante oblique orientée par la droite E = (Ipk = 1). Munis de ces notations, chaque composante de la distribution aléatoire s'exprime comme : fBetea,b)(X) r(a)r(b) x La distribution aléatoire f , y) s'exprimera alors : 3034221 18 f 9 k(E , E, y) = .W2(É> Ifik,Ek) sN' (ylig - g k (É>), e- el 7 e..- k(É)) où, .7«. Iµ, 62) représente la loi gaussienne de kt et de variance 62, Me Irt, E) la loi gaussienne bivariée de moyenne fi E rz2 et de matrice de covariance E. Une distribution gaussienne multivariée .7V'p(161/) est définie pour 5 X E 1 e-1(x - ii)' -Sr' - (x- g) fjv- oisi)(X) = 2 (27-0 p /21n11 /2 où IAI représente ici le déterminant de la matrice A. Note : cette définition couvre également le cas univarié quand p = 1. La mesure a priori Go du processus de Dirichlet est choisie de la façon suivante. La loi a priori du paramètre iik est une gaussienne, la 1 0 variance -ck est distribuée suivant une loi gamma-inverse, tpk suit une loi a priori de type Bernoulli et l'a priori pour les coefficients de régression coefficients /3k est une loi normale (gaussienne) multivariée de dimension IX(É)I. En utilisant la règle de Bayes, nous déduisons de l'expression de f (E, E', y) que : . i(yiE,E) = y Wk j'r2 (É> 1111<' lk ) .7V'( ylig - g k g), e-eTe.k(É)). It1I1=1W1 Neri 1, 11) \ 15 Finalement, en calculant l'espérance conditionnelle IE (yIE,É) à partir de l'expression de ffylE,E'), on obtient l'expression du modèle de réponse spectrale S(E,E') : , ,, .SlE,E ) = IE (yIE,É) = EX' .1 voovk j\r2(ÉI(Tt>k(-k) fl'k gk(É). Ei.iw1X2 E Yi, i Cette dernière expression permet de rendre compte de l'a priori 20 correspondant à l'approche bayésienne non paramétrique pour la réponse spectrale débruitée et interpolée.

3034221 19 Les caractères hétéroscédastiques et non gaussiens en tout (E, E') sont pris en compte dans le modèle de ffylE ,E') grâce au mélange de distributions gaussiennes .7^T (yliejl<' - gel e-PP'4(É)). A partir de ce choix de modèle probabiliste et des données 5 observées par simulation physique (Ei,E1i,yii), le procédé va estimer la loi a posteriori f (E, E',y1E1, E --- En, Ein,, y',) et l'espérance conditionnelle de cette loi a posteriori : .§(E,E') = IE yE,E',E1,E --- En, Ein,, représentant l'intensité dans le spectre d'une énergie de sortie E', pour toute énergie d'entrée E.

10 Pour le calcul exact de la distribution a posteriori le procédé utilise, par exemple, un schéma d'approximation Monte-Carlo par chaîne de Markov (MCMC) de type échantillonneur de Gibbs afin de générer des échantillons de la loi cible. Afin de permettre la procédure d'échantillonnage MCMC pour des objets de dimension infinie (infinité de composante du DP), le procédé 15 adopte une approche dite d'échantillonnage par tranche d'après Kalli et al. (Kalli M, Griffin J E and Walker SG 2011 Slice Sampling Mixture Models, Statistics and Computing, 21, 93-105), où seulement un nombre aléatoire fini K de composantes est nécessaire à chaque itération. Pour toute itération (t) de la procédure MCMC, il est ainsi possible 20 de générer une réponse spectrale débruitée S(E, Er. Pour T générations, la distribution a posteriori de la réponse spectrale est approximée par l'ensemble des tirages S(E,Er pour t = 1, ...,T, et la réponse estimée (moyenne a posteriori) s'exprime : T 1 S(E,E)TS(E , E)(t) t=1 Il est également possible de calculer de façon analogue l'écart-type a 25 posteriori et les intervalles crédibles à partir de la collection IS(E,Erl. Un exemple pour exécuter la procédure d'échantillonnage selon Gibbs va maintenant être donné.

3034221 20 Le modèle génératif pour les paramètres du mélange DPM est donné par : V1, V2, ... Beta(1, a) soient w1 = V1,w2 = V2(1 - V1), ... 111,112, --- J\r2(rTiii, -1,11, Gamma(ar, br) 1-1 1-1 T 1 T 2 Gamma(cy,b1-) --- Nix(É)1(m,g,r,g) 11)1,1)2, --- Bernoulli(p) soient Oi = il, pl, th), 02 = (µ2,i2, P2, 1P2), avec /1.1 = Tl = (Tl, ... Ce modèle génératif est entièrement équivalent à la densité de probabilité a priori f(wl,w2, 02, ...), et correspond à un a priori de type processus de Dirichlet pour la mesure aléatoire G(') = Ek's=i Wk 89k(') ^' D P(a, Go).

5 La description de la modélisation des paramètres DPM est complétée par un modèle génératif des données observées issues de la simulation physique (Ei,ri,yii). Pour cela, on introduit des variables latentes dites de classification KL1, définies pour i = 1, ...,n et j = 1, ...,n', tel que Ki = k si (Ei,ei,yii) est distribué suivant la kième composante du mélange 10 f(E,E,y). Boucle sur tous les n énergies d'entrée E et n' énergies E' de la grille de sortie Pour tout i < n,j < Kij I W1, W2, ... ^ Wk 81<e) k=1 01, - R7h p 1 1-Kij n avec r-Kij PKij.

0 T'Ici] 01, 02, (yiiI 161('ii - gKii(Ei, e PT (ii.5?1( ii(E j)).

3034221 21 Ce modèle génératif est entièrement équivalent à la vraisemblance : f (K11, , Knw, E1, , En, en,, v 11, --- v v1, 14/2, , 611, 02,.. A partir de ces distributions de probabilités et par application de la règle de Bayes, on cherche à calculer la densité de probabilité conditionnelle w2, 01, 02, ... I K11, , Knw, E1, , En, en:, Y11' --- Ynn'- 5 qui représente la distribution des pondérations et paramètres des composantes conditionnellement aux observations et aux variables de classification. Cette inférence est réalisée au moyen d'un échantillonneur de Gibbs qui à chaque itération (t) de l'algorithme, génère successivement les échantillons 10 suivants, pour tout k, 1411,1412, --- I K11, --- Knn' IPki ltk, 11-k, K11, --- Knw, , En, en, T kl K11, --- Knw, , En, en, K11, - - - Knw, , En, en, 12k' Tk,Tik,1Pk, K11, --- Knw, , En, en, (3 xi lek, ktik,Tk,Tik,1Pk, K11, --- Knw, , En, en:, Y11, - - - ,Ynn' et au niveau des nombres de points i dans la grille d'entrée d'énergie et de points j dans la grille de sortie, pour tout i < n, j < , PYij, [t2, --- il: vil: 2, --- T1, T --- T'2, --- ,I(31,1(32, --- ,1P1,1P2, --- W1, W2, --- Le comportement non paramétrique du procédé est caractérisé par la collection potentiellement infinie de paramètres wk et Bk. Afin de 15 réduire la difficulté calculatoire induite par le nombre infini de paramètres, il est possible d'utiliser une troncature du DP à un nombre suffisamment élevé de composantes (Ishwaran and James 2001). Une solution alternative proposée par (Kalli et al. 2011) est d'introduire une variable auxiliaire qui permet de ne générer qu'un nombre 20 aléatoire fini K de composantes à l'itération (t) tout en évitant une troncature arbitraire du modèle. Dans ce cas, les différentes étapes de génération aléatoire de l'algorithme seront les suivantes : 3034221 22 Initialisation aléatoire. - Générer pour 1 n',ui1-Uniform(0, 1/n . ni) , et poser u* = minu({util) - Générer V1-Beta(1, a) ; Poser w1 = V1 and r1 = 1 - V1. 5 - Pour k > 1, générer Vk-Beta(1, a), poser wk = Vkrk_i, poser rk = rk_1(1 - Vk), tant que rk > u*. - Affecter à K* la valeur maximale de k. - Générer 61k pour 1 < k < K*, à partir des distributions a priori. À l'itération t = 1, T, 10 - Générer (Ki./ 1E ii2, --- v 11'2, --- T1, T2, --- T'1, T'2, --- ,1(31,1(32, --- ,1P1,1P2, --- VV1, VV2, --- ) pour i n, j n' Calculer pour 1 < k < K*, 4 = (wk > ui1) max(wk, 1/n .0 Ek) g k , k( où 15 1(A) est la fonction indicatrice : 1(A) = 1 si A est vrai et 1(A) = 0 sinon. 1 re x Générer K - - - ric,*-12-ic k=1 - Réordonner les étiquettes de composantes en suivant leur ordre d'apparition dans la génération des Kij. Affecter à Kn le nombre de 20 valeurs distinctes de Kij. Poser, pour tout k Kn, nk = #{1(ii = k}, le nombre d'observations yii assigné à la composante k. - Générer (w1, w2 - - - , u11, - - - , unn' I K11, - - - , Knn') - Générer (W1, W2, wkn,rkn) Dirichlet(nl, n2, , nkn, a) - Générer pour i n', ut] -Uniform(0,min(wk, 1/n. n')) 25 - Poser u* = minu({util) - Pour k > Kn, générer Vk-Beta(1, a) , poser Wk = Vk rk_i, poser rk = rk_1(1 - Vk), tant que rk > u*. - Affecter à K* la valeur maximale de k. - Générer (lbI ki K11, --- Knn:, , En, E'n,) pour k , 3034221 23 bat. bai' (hr-FEr )(Er/202 t,j: Kii=k )(E'rifk)2 2 {ij: Kii=k )a -r-E( - Calculer vo = nk nk - Calculer = btat. bt,at' nk nk )a.c+ )cit-F( 2 2 (bt-*{i,j: +(erifk) (Erilk)(erik) bt41E{ Kii=k)(Ei-iik)2+(erg'k) +2(E i- iik)(E'r g'k) ti - 5 - Générer v- Uniform(0,1) - Si < , poser tpk = 0 sinon poser tpk = 1 vo+vi - Générer (-cklk,---, Knw, , En, en') pour k < K*, - Si 11)k = 0 nk -c,7 1-Gamma + -2, + -12 (E -11x)2 - Si th< = 1 T71-Gamma I c, +7nk , br + 4 1 (Ei - µk)2 + (ri - µ'k)2 - 2 (Ei - uk)(Ei - µ'k)\ { i,j: Kii=k } i 10 - Générer (Tikl ktk, [t k, lPk, K11, --- , Knn', El, el, --- , En, en') - Si tpk = 0 T -Gamma 2 \ nk 1 v, +, br, + -2 (ri - - Si th< = 1 / \ -1-'71-Gamma ar- ++ 'hi- + il 1 (Ei - ilk)2 + (ri - Il' k)2 + 2 (Ei - PO (ri - iek) fi,j: Kii=k1 i - Générer (fikl T lb -k, - T1k, T k, - K 11, --- , Knn:, Ei, el, ... , En, E n') pour k < K*, T - Calculer Sitic = (E{ii: K ii.k} E i' Etii:Kii=k1E;) 15 - Poser Ek = Rip- kl.(ok TO,). Rek avec Rek = (0 11 °) si tpk = 0 et Rek = .'7(11 -11) si lPk = 1 =k 1 pour k < K*, 3034221 24 - Calculer Vµk = nkEk-1)-1 - Calculer Mil* = Vitk.(1-1:1 - Générer µk-s7V'2(Mitk, Viik) - Générer I (8 k. Lt k,Tk,Tk,K,K11,---, Knn,, --- , En, E yll, --- , ynn,) 5 pour k < K*, - Poser rflk fek T Kii=k1 - Xk - Poser fek( Kii=k1 Mflic = rflk - - yii - e-Yii ri-gl - Mig - Générer lek - .Np-e(É)1 (Mflk, rigk). Pour ce tirage aléatoire, on approxime la variance d'observation 10 de yii (elle-même approximée à une variable aléatoire gaussienne) = e-flii-4(Éii) à la valeur de l'observation Poisson Âij = e-Yii correspondante. Cette approximation, qui simplifie grandement l'inférence, est d'autant plus valable que Âii est élevé. - Calculer à chaque itération (t) pour tout choix de grille de prédiction 15 É = (E,E), y 1E S (E, Eo =(yIE ,E") Wk Neri k,zk) k k=iri=iwi pzi) le - xk(E) La réponse spectrale débruitée peut, au choix, partager ou non la même grille (Ei, ri) que la simulation Monte-Carlo physique initiale. En effet, l'utilisateur peut choisir n'importe quel point (E, E) d'intérêt. En particulier, il est possible d'interpoler entre les points de la grille initiale.

20 Finalement, à partir de T itérations de l'échantillonneur de Gibbs, on obtient une estimation du méta-modèle pour la spectrométrie 3034221 25 T (E,E1) -11S (E,E1)(t) t=i Il est alors possible de calculer l'écart-type a posteriori de l'estimateur de réponse spectrale G (S (E,E' ) (E, E') - S (E, 012) .7 T -1 1 Les hyper paramètres peuvent être considérés fixes dans l'estimation de la réponse du système spectrométrique en fonction de 5 connaissances physiques préalables sur le détecteur. D'un autre côté, ceux- ci peuvent être estimés au moyen d'un niveau de hiérarchie supplémentaire qui permet de leur attribuer un a priori dit vague. - À titre d'exemple, il est possible de placer un a priori Gamma@Pb, sur les paramètres d'échelle br et b", dans la distribution des amplitudes de 10 composantes. Ce choix conduit à une distribution a posteriori Gamma. Exemple : br- Gamma Kn 1 (Pb + Kn.,41 +- k-Ck =1 - Il peut être pertinent d'estimer le paramètre de concentration a du processus de Dirichlet en suivant l'approche de Escobar and West (1995). 15 - Le paramètre p de la loi de Bernoulli intervenant dans le modèle de mélange du type de composante tpk peut être lui-même estimé par l'intermédiaire d'un hyperprior Beta(rt-0,71-1). La loi a posteriori pour p suit alors l'expression : p- Beta(zo + #{k < Kn 11)k = 0},71-1+ #{k Kn: 11)k = 1}) - Notons également qu'il peut être efficace de remplacer l'a priori gaussien 20 sur les coefficients de régression /3k par un a priori dit empirique basé sur le tirage aléatoire de P points (É1,E,_,...,Ép,E"p) à partir de .7V'2(fik,Ek) et de prendre pour (5,1, la plus proche valeur de yii correspondant à t=1 3034221 26 chaque point échantillonné. P doit être plus grand que la taille du vecteur pk. Nous générons alors lek - .7V'(Mig, rg) comme loi a priori (avec Ép = (4, Éip )) : P 1 \- rek T . e- 571) rig = - X \p=1 Mfl = rig - (Ei)c=i Xkgp) Sip e-51'19), 5 L'expression pour la loi a posteriori de lek en découle directement. A l'issue de ces différentes étapes, le procédé dispose d'un méta-modèle continu caractérisant tout type de détecteur pour la spectrométrie nucléaire, le méta-modèle prenant en compte l'ensemble des interactions physiques intervenant dans la réponse spectrale pour une énergie d'entrée 10 donnée. Le méta-modèle ainsi généré pourra être injecté dans une méthode de convolution de spectre nucléaire utilisant une modélisation probabiliste non paramétrique telle que celle proposée par Barat et al 2007 (Barat, E.; Dautremer, T.; Montagu, T., "Nonparametric bayesian inference in 15 nuclear spectrometry," Nuclear Science Symposium Conference Record, 2007. NSS '07. IEEE , vol.1, no., pp.880-887, Oct. 26 2007-Nov. 3 2007), afin de déterminer les radioéléments présents dans le spectre et leur activité. En fonction de l'environnement dans lequel se trouve l'objet à analyser, des variations de forme de signature spectrale dues à la présence 20 d'effets de matrice au niveau de la source d'entrée peuvent apparaître. Le méta-modèle défini selon l'invention peut prendre en compte plusieurs réponses spectrales estimées comme il a été explicité ci-avant. Pour cela, une solution consiste par exemple à utiliser différents blindages de différentes épaisseurs et de différents matériaux, pour être simulés et 25 intervenir dans le méta-modèle. Une autre solution consiste à conserver un arbre de Pàlya pour modéliser les interactions diffuses conduisant à des fonds continus qui pourraient ne pas être modélisés dans la réponse.

3034221 27 Dans l'hypothèse où le méta-modèle prend en compte une collection d'effets de matrice comme précédemment décrit, le spectre d'entrée recherché n'est alors constitué que de pics discrets. Il est alors possible de substituer à l'a priori par processus de Dirichlet sur des raies 5 mono-énergétiques, un a priori par processus de Dirichlet sur les radioéléments d'une base de données de radionucléides. Le spectre recherché est alors constitué d'une somme discrète d'éléments du tableau périodique dont le nombre de composantes est illimité. Il devient alors possible d'affecter à chaque radioélément une probabilité a priori de 10 présence dans l'échantillon analysé. Le méta-modèle est à nouveau utilisé pour caractériser la réponse spectrale du détecteur de l'élément considéré. Cette réponse est directement déterminée par la somme discrète des réponses spectrales individuelles correspondant à chaque raie mono-énergétique considérée. L'utilisation d'une base de données de 15 radionucléides permet notamment de considérer la calibration en énergie du spectre observé comme incertaine et d'estimer celle-ci à partir des données. Pour chaque valeur proposée de keV/canal il est nécessaire d'évaluer la réponse spectrale du système pour chaque énergie de la grille de sortie. Cette opération utilise le méta-modèle continu, éliminant ainsi tout besoin 20 d'interpolations dans le cas d'un recours à un modèle discret. Dans le cas où le méta-modèle prend en compte une collection d'effets de matrice, pour estimer directement l'activité des radionucléides, on utilise un méta-modèle étendu noté S(E,E'', 0 où représente un paramètre caractéristique de l'effet de matrice associé à la source. La collection d'effets 25 de matrice est discrète de taille M et chaque élément est identifié par un indice entier. L'algorithme d'estimation des activités est fondé sur une approche Monte Carlo par chaîne de Markov (MCMC) et plus particulièrement sur un échantillonneur de Gibbs dont un exemple est détaillé ci-après.

30 Avant d'expliciter la description de l'algorithme, différents paramètres sont posés. Le spectre est considéré comme un mélange par 3034221 28 processus de Dirichlet de différentes réponses de radionucléides normalisées. Un exemple de construction de ce mélange est donné ci-après. On note x un radionucléide quelconque d'une table donnée de radioéléments. Cette table étant discrète, on repèrera dans la suite chaque 5 élément par un indice entier. On note NX le nombre d'émetteurs retenus pour l'élément x considéré et ifx,/ pour / = 1, ...,Nx, les probabilités d'émissions associées ainsi que vx,/ pour / = 1, ...,NX, les énergies correspondantes. On définit la réponse du radionucléide tpx(E, 0, pour une énergie observée E' et un effet de matrice par : NX tpx,-(E') S(v xj, E, 0, 1-1 * Tcx,i 10 on définit = fo tpx,e1) dE' et rcx/ = pour tout 1 = 1, , Nx. On définit alors, une réponse de radionucléide normalisée par : NX tpx* ,(E1) =n-x* 1=1 on vérifie que fo- tpx* e1) dE' = 1. La densité de probabilité du ième photon observé dans le canal d'énergie E; s'exprime alors pour i = 1, ...,n: J (Ei) = Wk E;) k=1 15 où p est un paramètre positif de conversion de l'indice de canal en énergie (keV/canal). Nous suggérons pour p un a priori gaussien centré sur ktp et d'écart-type a : P J'r (P I Pp, Ce mélange par processus de Dirichlet s'appuie sur la mesure de probabilité F générée suivant un processus de Dirichlet : F DP(a, FXFÎ) 3034221 29 FX (x) = 8i(x) représente la distribution de probabilité a priori d'observer un élément x où J est le nombre de radionucléides de la base retenue pour l'analyse et pf > O. Fî est la loi a priori uniforme discrète sur la collection d'effets de matrice, et a est le paramètre (positif) de concentration 5 du processus : F(-) = Wk 8X/Jk(*) k=1 avec w1 = V1, wk = Vk fi i (1- V1) tel que Vk-Beta(1, a) et (Su(-) représente la fonction de Dirac localisée en u. Avant de décrire les différents tirages aléatoires permettant d'explorer la loi a posteriori, il reste à introduire les variables Ki d'allocation du ième photon 10 observé à une composante k du mélange. L'échantillonneur de Gibbs consiste à alterner des tirages aléatoires suivant les lois conditionnelles suivantes, pour tout i < n, w1,w2, --- X1, X2, --- ---, P et pour tout k 1 Ki, Kn Xk I Kn, E1, - - - P kI K1, --- Kn,E;, , En' , Xi, X2, P ainsi que pl K1, , Kn, , En, X1, X2, ---, --- 15 Comme pour la détermination du méta-modèle, afin de minimiser le nombre de calculs, l'approche proposée par (Kalli et al. 2011) et l'introduction d'une variable auxiliaire ui qui permet de ne générer qu'un nombre aléatoire fini K de composantes à l'itération (t) tout en évitant une troncature arbitraire du modèle sont utilisées. Sous ces hypothèses, les 20 différentes étapes de génération aléatoire de l'algorithme sont détaillées ci- après. lnitialisation aléatoire. - Générer pour 1 ui-Uniform(0, lin), et poser u* = minu({uil) 3034221 30 - Générer V1-Beta(1, a) ; poser w1 = V1 et r1 = 1 - V1. - Pour k > 1, générer Vk-Beta(1, a), poser wk = Vk rk_i, et poser rk = rk_1(1 - Vk), tant que rk > u*. - Affecter à K* la valeur maximale de k. 5 - Générer Xk pour 1 < k < K*, à partir de la distribution a priori FI). - Générer k pour 1 < k < K*, à partir de la distribution a priori - Générer p (plup, À l'itération t = 1, T, - Générer (KilE;, w1, w2, X2, p) pour i < n, 10 - Calculer pour 1 < k < K*, - = l(Wk > ni) MaX(Wk, IP x* k,fic(PE;) où 1(A) est la fonction indicatrice : 1(A) = 1 si A est vrai et 1(A) = 0 sinon. - Générer Ki- -El<-1 8k(Ki), - Réordonner les étiquettes de composantes en suivant leur ordre 15 d'apparition dans la génération des Ki. Affecter à Kn le nombre de valeurs distinctes de Ki. Poser, pour tout k Kn, nk = #{Ki = k}, le nombre d'observations E; assignées à la composante k. - Générer (wi, w2 ..., uni Kn) - Générer (wi, w2, , wK',r,') Dirichlet(ni, n2, , nicn, a), 20 - Générer pour i n, ui-Uniform(0 ,min(wk, 1/n)), - Poser u* = minu({uil), - Pour k > Kn, générer Vk-Beta(1, a), poser Wk = V r Vk k-15 poser rk = rk_1(1 - Vk), tant que rk > u*, - Affecter à K* la valeur maximale de k. 25 - Générer (xk I Kn, p) pour k < K*, - Calculer pour tout j = 1, ...,J, = i infic(PE;), Ki=k 1 - Générer xk- 1 EJ, Mxk). Ei=ink.l .k,j 3034221 31 - Générer (i<1 ...,Kii,E;, ..., ,x1, x2, , p) pour k , Calculer pour tout m = 1, ...,M, = tpx*k,m(pEi), K i=k Générer Lin.1 ri2.1 1<,m 8mgk). - Générer Coi K1, , Kn, Ev...,E Xl, X2, - - - - - -) par une étape 5 Metropolis-Hastings, Générer une marche aléatoire d'écart-type Ep autour de p P* J'r (P* IP, En. - Calculer N(P*111^0' nriL=1IPXKi,Ki(CP*E;) eMH J'r (PIPP' 1 IP,(KI,Ki(PEi) Générer une variable uniforme y - Uniform(0,1).

10 Si y min (1, emH) affecter p = p*, sinon laisser p inchangé. Finalement, à partir de T itérations de l'échantillonneur de Gibbs, on obtient une estimation des activités des radioéléments intervenant dans le mélange, pour tout k, T A- 1 wl,(,t)A(t) Xk T - t=1 De surcroît, il est également immédiat de calculer l'écart-type a posteriori des 15 activités : 1 T AXk (T - 1 Li \ Xk k 1 Gr w(t)A(t), )2 Y Il est également possible d'obtenir le spectre d'entrée estimé, déconvolué de la réponse du système N (t) T Xk (t) S(E) Wk IrXk (t) 1 (t) (5v (t) (E) Xk '1 t=1 1=1 t=1 )2 3034221 32 Dans la version proposée l'algorithme permet l'estimation du gain de conversion keV/canal à partir des seules données ainsi que la détermination d'un effet de matrice (potentiellement différent) pour chaque élément intervenant dans le mélange. Il permet également de fixer une 5 probabilité a priori d'occurrence des radionucléides dans le contexte considéré et d'obtenir directement les activités des éléments constituants le mélange ainsi que les incertitudes associées. Les figures 9 et 10 illustrent le spectre des énergies obtenu à partir de l'analyse du spectre de la figure 10 en utilisant le procédé selon 10 l'invention. La courbe en pointillé représente le spectre simulé, la courbe en trait plein le spectre reconstruit. Les figures 11 et 12 illustrent un deuxième exemple de spectre. La figure 11 représente le spectre des énergies incidentes, avec une région d'émission du premier pic au 60Co, 1173,2 keV. La figure 12 le même spectre 15 avec un zoom sur la région d'intérêt.

Claims (13)

1. REVENDICATIONS1 - Procédé pour déterminer la nature des radioéléments présents dans un objet et leur activité caractérisé en ce qu'il comporte au moins les étapes suivantes : - une première phase (32) de simulation numérique de réponses spectrométriques pour un ensemble d'énergie incidente E et un ensemble d'énergie de sortie mesurée E', afin d'obtenir un ensemble de données simulées, - une deuxième phase de régression non paramétrique (33) sur les données simulées, estimation non paramétrique de la quantité représentant la probabilité jointe des triplets (E,E',y) à partir de points simulés (Ei,Ei',yij) afin d'en déduire un méta-modèle S(E, E') pour tout couple d'énergie (E, E') sur une fonction continue, - à partir du méta-modèle S(E, E'), (36) la détermination de la nature et de l'activité des radioéléments présents dans l'objet.
2. 2 - Procédé selon la revendication 1 caractérisé en ce que : - on utilise n un nombre de points dans la grille d'entrée, énergies E, et n' un nombre de points dans la grille de sortie, énergies E', - pour i = 1, ...,n et j = 1, ...,n', on dispose des données caractéristiques des intensités spectrales calculées À pour une énergie d'entrée EL et une énergie de sortie E'1, - on utilise une méthode d'estimation non paramétrique de la quantité f(E,E',y), représentant la densité de probabilité jointe des triplets (E,E',y) à partir des points simulés (Ei,E'pyii), et on déduit un modèle S(E, E') pour tout (E, E') E où I': est un espace continu : rz 3034221 34 S(E,É) = IE (yIE,É) fR y - f , , y) dy = I y - f (yIE , dy R fR f , y) dy où le symbole lE(y) représente l'espérance mathématique de la variable aléatoire y, y est déduite des intensités spectrales calculées xv. 5
3. 3 - Procédé selon la revendication 2 caractérisé en ce que l'on détermine les valeurs À en utilisant un logiciel Monte-Carlo, les données À étant considérées comme des réalisations issues d'une distribution de Poisson dont l'intensité ILS est estimée au moyen d'une procédure de régression non paramétrique et on introduit yii = + E) où 0 < E « 1, puis 10 on approxime la distribution de probabilité de yii pour des valeurs suffisamment grandes de IL> > 10 par une loi gaussienne de moyenne log ILS et de variance L pour la loi jointe f , , y) , on choisit un mélange par processus de Dirichlet (DPM) comme distribution a priori et on exprime la distribution aléatoire en 15 une somme sur l'infini de foi, composantes de f , , y), f(E, y) =wk fok(E, y) k=i paramétrée par 61k le paramètre associé à la kième composante de G mesure aléatoire définie par G(.) = Wk 8ok k=1 avec w1 = V1, Wk = - V1) tel que Vk -Beta(1, a) et (Su(-) représente la fonction de Dirac localisée en u et Beta(a, b), pour 0 < x < 1 avec r(a+b) 20 fBeta(a,b) (x) = nor(ox a-1 b-i (1 - x) . 3034221 35
4. 4 - Procédé selon la revendication 3 caractérisé en ce que les composantes de la loi jointe foi, sont exprimées à partir des valeurs suivantes : - 0k = T k,ipk,f3k) avec fik =µlkl = lek, k), - Xk(E) le vecteur de régresseurs centré avec É = (E,E'), et /3k le vecteur 5 de coefficients de régression, - la matrice Ek, dépendant du paramètre tpk E {0,1}, comme suit : Ek = ek 2 1 = (1 1) si R-1 (Tk ° R avec R 1- ( °) si tpk = 0 et R ek. -ci). 0 1 lPk = 1, permettant de choisir entre une composante alignée sur les axes E et E' (Ipk= 0) et une composante oblique orientée par la droite E = 10 (Pk = 1). chaque composante foi, s'exprime alors : fek(E, E , y) = .7V'2(É> I fik, k) N (ylig - g k(É>), 0-2) représente la loi gaussienne de kt et de variance 62, .N2(- E) la loi gaussienne bivariée de moyenne fi E rz2 et de matrice de covariance E 15 où la loi a priori du paramètre µk est une gaussienne, la variance -ck est distribuée suivant une loi gamma-inverse, tpk suit une loi a priori de type Bernoulli et l'a priori pour les coefficients de régression coefficients /3k est une loi normale (gaussienne) multivariée de dimension IX(É)I, et 20 en appliquant une règle de Bayes sur l'expression de ffE,E',y) on obtient l'expression f fylE , E = w k j'r2 (É> fi 11:1' k) \ (ylig - g kg), e- el( k(É)) J2 (E irt1,zo et le modèle probabiliste S(E ,E') = E (yIE,E') wk (É> 111:i' \ le' - g (É) wi J\r2(r El) k k 3034221 36 A partir de ce modèle probabiliste et des données observées par simulation Épyii) , on estime la loi a posteriori f(E,e,y1E1_,E Y11 --- En, Ei Ynni) et l'espérance conditionnelle .§"(E, E') = IE (y1E,E',E1_,E v --- En, Ei Ynni) afin de déterminer les éléments présents dans l'objet et leur activité. 5
5. 5 - Procédé selon la revendication 4 caractérisé en ce que pour le calcul de l'a posteriori, on utilise un schéma d'approximation Monte-Carlo par chaîne de Markov (MCMC), pour toute itération (t) de la procédure MCMC, on génère une réponse 10 spectrale débruitée S(E,E ft), pour T générations, la distribution a posteriori de la réponse spectrale est approximée par l'ensemble des tirages S(E,É)(t) pour t = 1, ...,T, et la réponse estimée s'exprime : T 1 T 1S(EEft) t=1 15
6. 6 - Procédé selon la revendication 5 caractérisé en ce que le schéma d'approximation comporte une étape d'échantillonnage par tranche utilisant un nombre aléatoire fini K de composantes pour chaque itération et en ce qu'il comporte les étapes suivantes : on introduit des variables latentes de classification Kii, définies pour i = 20 1,...,n et j = tel que Ki= k si (Ei,ri,yii) est distribué suivant la kième composante du mélange f(E,E",y), on définit un modèle pour les paramètres du mélange, pour tout i <n,j < KijIW1,W2, -- - 8k (-) k=1 01, 02, N2(Ei,E1i1,11Kii,Exi1) 3034221 37 avec LI p-i ° Ki] - 0 Ti - Ki] Kij, EL, 01, 02, (Yiji -XKiJ(EL,ej), e ifj?Kii(Ei"E'l)) équivalent à la vraisemblance f(K11, Knn', E1, En, en', Y11, - - - , Ynn'i w1, w2, - - - , 01, 02, - - -) à partir de ces distributions de probabilités et par application de la règle de Bayes, on calcule la densité de probabilité conditionnelle fW1, W2, - -- 611, 02, --- I K11, --- Knn', - - - , En, en', Y11, - - - , Ynn') 5 en utilisant un échantillonneur de Gibbs qui à chaque itération (t) génère successivement les échantillons suivants, pour tout k, w1, w2, --- I K11, - -- Knn' kI -- - , Knn', E - -- , En, E ktk, -- - , Ei,-- - , En, E ktk, [t k, Pk, K11,-- - , -- - , En, en' 11-1 ki k, k, k, K11, -- - , Knn', el, - - - , En, en' (3 xi 11-1 k, Tk,Tik,IN, K11, - -- Knn',E1,Ei,- - - , En, E in', Y11, - - - , Ynn' et au niveau des nombres de points i dans la grille d'entrée d'énergie et de points j dans la grille de sortie, pour tout i < n, j < E'1, Y ij, kt2, -" 1, 11:2, --- T1, T2, --- T '1, T '2, --- ,1(31,1(32, --- ,1P1,1P2, --- VV1, VV2, --- à partir de T itérations de l'échantillonneur de Gibbs, on obtient une 10 estimation du méta-modèle pour la spectrométrie : 1 (E, Ei) (E,E1)(t). t=1
7. 7 - Procédé selon la revendication 6 caractérisé en ce que l'on introduit une variable auxiliaire uij afin de ne générer qu'un nombre aléatoire fini K de composantes à l'itération (t) tout en évitant une troncature arbitraire du 15 modèle. 3034221 38
8. 8 - Procédé selon l'une des revendications 5 à 7 caractérisé en ce qu'il comporte une étape de calcul de l'écart-type a posteriori et des intervalles crédibles à partir de l'ensemble {S(E,E')(t)} 1 7 65 ,E) 1 , E') - S , E')(t))2) . 5
9. 9 - Procédé selon l'une des revendications 2 à 8 caractérisé en ce que l'on introduit un a priori Gamma4b, sur les paramètres d'échelle br et b'T dans la distribution des amplitudes de composantes conduisant à une distribution a posteriori Gamma : br- Gamma / Kn avec Gamma(a, b), pour x > 0, 1 (Pb Kn, +- k-Ck =1 1, b a xa-1 e-bx. fGamma(a,b) (X) r(a) 10
10. 10 - Procédé selon la revendication 4 caractérisé en ce que l'on remplace l'a priori gaussien sur les coefficients de régression lek par un a priori basé sur le tirage aléatoire de P points (Ê''1, à partir de .7V'2(fik, Ek) et on prend pour (5,1, la plus proche valeur de yii correspondant à chaque 15 point échantillonné où P est plus grand que la taille du vecteur pk, on génère alors /3k - .7V'(Mig,Fig) comme loi a priori (avec Ép p)) / P \1 r= g k gk(Ep) - e -5719 \p=1 me = re - L gk( \P.1 t=1 3034221 39
11. 11 - Procédé selon l'une des revendications précédentes caractérisé en ce que l'on génère un méta-modèle étendu noté S(E,E'',) où est un paramètre identifié par un indice entier et caractéristique d'un effet de matrice. 5
12. 12 - Procédé selon l'une des revendications précédentes caractérisé en ce que l'on estime l'activité des radioéléments en exécutant les étapes suivantes : soit Nx le nombre d'émetteurs retenus pour l'élément x considéré et n-x,/ pour 10 1 = 1, ...,Nx, les probabilités d'émissions associées ainsi que vx,/ pour 1 = 1, ...,Nx, les énergies correspondantes, on définit la réponse du radionucléide tpx(E, 0, pour une énergie observée E' et un effet de matrice par Nx tp =71" S(v xj, E', 1=1 * X.on définit = fo tp xe') de et rcx/ = - pour tout 1 = 1, ...,Nx, x,e 15 on définit une réponse de radionucléide normalisée par Nx tpx =71-x* S(v xj, E, 1=1 on vérifie que fo- tpx* (El) dE' = 1 on exprime la densité de probabilité du ième photon observé dans le canal d'énergie E; pour i = 1, n J (Ei) = Wk E;) k=1 où p est un paramètre positif de conversion de l'indice de canal en énergie 20 (keV/canal). On introduit les variables Ki d'allocation du ième photon observé à une composante k du mélange, 3034221 40 On alterne des tirages aléatoires suivant des lois conditionnelles suivantes, pour tout i < n, pour tout k wi, w 2, ... 1 Ki, ... , Kn hl K1, ... , Kn, E1, ... , En, 1, .2, ---, /0 .kI K1, ... ,Kn, E;, ... , En' , Xi, X2, ... , P ainsi que pl K1, ... , Kn, E,_, ... , En , X1, Y 2, ---, 1, 2, --- 5 en utilisant un nombre aléatoire fini de composantes à l'itération (t), à partir de T itérations de l'échantillonneur de Gibbs on obtient une estimation des activités radioéléments intervenant dans le mélange pour tout k, T A- 1 ',k,' w(t)A(t) A(t) Xk T k Xk,-,,,- k t=1 avec p un a priori gaussien centré sur iip et d'écart-type ap. P - -WC° I Pp, c)-i, ).
13. 13 - Procédé selon la revendication 12 caractérisé en ce que l'on calcule l'écart-type a posteriori des activités GAXk ''' T w(t)A(t) )2 1 k Xic,k )2 ( 1 Y ('4--- V - 1 L \ Xk _i t=1 et/ou le spectre d'entrée estimé, déconvolué de la réponse du système N (t) T Xk 1 ( S(E) R.--,' T Wkt) ITXk 1 (t) (t) (5v (t) (E) . ,,- /c Xk 't t=1 1=1 10 15