EP3274863A1 - 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

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 ¹³⁷Cs émettant un gamma à 661 ,5 keV ou le ⁶⁰Co émettant un gamma à 1 173,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 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(l) obtenu à l'aide d'un détecteur HPGe et d'une source de ¹³⁷Cs é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.

II 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 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(ll) obtenu à l'aide d'un scintillateur plastique de type EJ200, en présence de sources ⁶⁰Co et ¹³⁷Cs. Les formes caractéristiques de ces spectres sont dues aux fronts Compton 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 ¹³⁷Cs et de ⁶⁰Co. 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. 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 radioéléments émetteurs et mesurer l'activité de ces derniers, consiste à analyser le spectre dans son ensemble et à ne pas se restreindre 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 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 é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 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 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 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 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 simulation de la réponse d'un scintillateur plastique de type EJ200, en présence de trois sources gamma (²⁴¹Am, ¹³⁷Cs, ⁶⁰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 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 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 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.

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. 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 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 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 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 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 ⁶⁰Co on simulera directement la réponse du détecteur pour deux énergies incidentes, émises respectivement à 1 173,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 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.

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 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 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 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'), 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 n' un nombre de points dans la grille de sortie, énergies £",

• pour i = l, ... , n et j = l, ... , n', on dispose des données caractéristiques des intensités spectrales calculées À_£y- pour une énergie d'entrée £_£ et une énergie de sortie E ,

• 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 (Ει, Ε , γ^), et on déduit un modèle S(E, £") pour tout (E, £") E KL² où IL est un espace continu :

où le symbole E(y) représente l'espérance mathématique de la variable aléatoire y, y est déduite des intensités spectrales calculées λ ·.

On peut déterminer les valeurs λ_£;- en utilisant un logiciel Monte-

Carlo, les données X étant considérées comme des réalisations issues d'une distribution de Poisson dont l'intensité est estimée au moyen d'une procédure de régression non paramétrique et on introduit y_tj = log(À_£y- + ε) où 0 < ε « 1, puis

on approxime la distribution de probabilité de y_tj pour des valeurs suffisamment grandes de > 10 par une loi gaussienne de moyenne log et de variance , pour la loi jointe f(E, E' ,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 f_Q]i composantes de f{E, E ,y), f(E, E^', y) = ^ w_k fe_k (E, E^',y)

k=l

paramétrée par 9_k le paramètre associé à la k^ième composante de G mesure aléatoire définie par

G(-) = w_fc <¾(·)

fc=l

avec w₁ = V_x, w_k = V_k a) et 5_U(-) représente la fonction de Dirac localisée en u et Beta(a, è), pour 0 < x < 1 avec Les composantes de la loi jointe f_Qk sont, par exemple, exprimées à partir des valeurs suivantes :

• 0k = (fik> fk>≠k k avec ¾ = (μ¾, μ'_¾), ¾ = (τ_¾,τ'_¾),

• ¾ (E) le vecteur de régresseurs centré avec È = (Ε, Ε^') , et /?_fe le vecteur de coefficients de régression,

· la matrice ∑_fe, dépendant du paramètre ifj_k G {0,1}, comme suit : ∑_k =

= f (1 ^"l¹) ⁵¹

= 1, permettant de choisir entre une composante alignée sur les axes E et E^' {ifj_k = 0) et une composante oblique orientée par la droite E = E Wk = !)^■ chaque com osante f_9k s'exprime alors :

où, Χ ·\μ,σ²) représente la loi gaussienne de μ et de variance σ², Κ₂{· \μ,∑) la loi gaussienne bivariée de moyenne μ G R² et de matrice de covariance∑, où

la loi a priori du paramètre _¾ est une gaussienne, la variance r_k est distribuée suivant une loi gamma-inverse, ifj_k suit une loi a priori de type Bernoulli et l'a priori pour les coefficients de régression coefficients _k est une loi normale (gaussienne) multivariée de dimension |^(E)|, et

en appliquant une règle de Bayes sur l'expression de f{E,E,y) on obtient l'expression

et le modèle probabiliste

à partir de ce modèle probabiliste et des données observées par simulation

(Ei.E'_j.yi_j) , on estime la loi a posteriori f E,E' ,y\E_x,E ,y_xx, ...,E_n,E'_n y_nni) et l'espérance conditionnelle S E,E') = E (y\E,E' ' ,E_X,E' ,y_xl, ...,E_n,E'_n,,y_nni) 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')^, pour T générations, la distribution a posteriori de la réponse spectrale est approximée par l'ensemble des tirages S{E,E')^ pour t = 1, ... , T, et la réponse estimée s'exprime :

t=i Le schéma d'approximation peut comporter une étape d'échantillonnage par tranche utilisant un nombre aléatoire fini κ de composantes pour chaque itération et en ce qu'il comporte les étapes suivantes :

on introduit des variables latentes de classification K^, définies pour ί = Ι,.,.,η et ; = Ι,.,.,η', tel que K_t = k si (E_£ est distribué suivant la fciè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 < n,

équivalent à la vraisemblance

f(Kll_>■■■ _> Knn'_> E_\, E _1} ... , E_n, E _n;

à partir de ces distributions de probabilités et par application de la règle de Bayes, on calcule la densité de probabilité conditionnelle

f(w₁,w₂, ...,θ , θ₂, ... I K_llt ... ,K_nn>,E_x,E , ...,E_n,E'_n>,y_xx, ...,y_nn en utilisant un échantillonner de Gibbs qui à chaque itération (t) génère successivement les échantillons suivants, pour tout k,

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≤ n,

ι,γιι_>μ₁,μ₂, ...,μ μ ₂, ...,τ₁,τ₂,—,τ ₂,—,βι,β2>■■■>Ψι>^'Ψ2>■■■ >^wi>^w2>■■ à partir de T itérations de l'échantillonneur de Gibbs, on obtient une estimation du méta-modèle pour la spectrométrie

T t=i

Il est possible d'introduire une variable auxiliaire afin de ne générer qu'un nombre aléatoire fini κ 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')^}

Selon une variante de réalisation, on introduit un a priori Gamma(<p_fc,¾) sur les paramètres d'échelle b_T et b'_T dans la distribution des amplitudes de composantes conduisant à une distribution a posteriori Gamma :

avec Gamma(a,ô), pour > 0,

b^a _a_₁ __bx

Gamma(a,fc) 0Ό^— |^~(a) ^ ^

L'a priori gaussien sur les coefficients de régression β_1ί peut être remplacé par un a priori basé sur le tirage aléatoire de P points (Ei,E'i, ...,Ë_P,E'p) à partir de ^\₂( _¾,Σ_¾) et on prend pour ( _1; ...,y_P) la plus proche valeur de y_tj correspondant à chaque point échantillonné où P est plus grand que la taille du vecteur /?_fe, on génère alors _k ~ Ν(Μ_β,Γ_β) comme loi a priori (avec Ë_p = (Ë_P,Ê'_P))

Selon une variante de réalisation, le procédé génère un méta- modèle étendu noté Ξ(Ε,Ε^',ξ) où ξ est un paramètre identifié par un indice 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 Ν_χ le nombre d'émetteurs retenus pour l'élément χ considéré et π_χ pour l = l,...,N_x, les probabilités d'émissions associées ainsi que v_x pour 1 = 1, ...,Ν_χ, les énergies correspondantes,

on définit la réponse du radionucléide -ψ_χ{Ε',ξ), pour une énergie observée E' et un effet de matrice ξ, par

1=1

on définit^ = J₀ Ψχ_,ξίΕ¹) dE'et π^* _ιξ = ^ pour tout / = 1,...,Ν_χ, on définit une réponse de radionucléide normalisée par

on vérifie que J₀ ¾(£") dE' = 1

on exprime la densité de probabilité du i^eme photon observé dans d'énergie pour i = Ι,.,.,η f{E[) = _jWk r_Xk,_?k{p-^Ei^')

k=l où p est un paramètre positif de conversion de l'indice de canal en énergie (keV/canal).

On introduit les variables K_t d'allocation du i^eme photon observé à une composante k du mélange,

On alterne des tirages aléatoires suivant des lois conditionnelles suivantes, pour tout i≤ n,

pour tout k

w_ltw₂, ... I K_lt ...,K_n

Xk \ Κ-L, ...,K_n,E_lt ... , Ε_η,ξ₁,ξ₂, ...,p

Κι_>—,K_n,E , -,Ε_η,χ₁,χ₂, ...,p

ainsi que

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,

avec p un a priori gaussien centré sur μ_ρ et d'écart-type σ_ρ.

On calcule par exemple l'écart-type a posteriori des activités en exécutant les étapes suivantes :

et/ou le spectre d'entrée estimé, déconvolué de la réponse du système 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 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,

• 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,

• 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 œuvre du procédé selon l'invention,

• La figure 8, un synoptique des étapes mises en œuvre par le procédé 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 1 1 et 12, un deuxième exemple de spectre des énergies incidentes obtenu par la mise en œuvre de l'invention.

La figure 7 est un exemple de système permettant la mise en œuvre 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 , relié à une électronique d'acquisition permettant d'obtenir un spectre, un processeur 23 ou unité informatique adapté à exécuter les étapes du 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 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 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 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.

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 exemple de la connaissance physique, 34, sous forme d'à 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 constantes quelle que soit l'énergie d'entrée E,. A partir du méta-modèle 35, et de la connaissance a priori, il sera possible de remonter, par application 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.

Le procédé va chercher à estimer une fonction continue de E² en

E et E' à partir d'un nombre potentiellement infini de paramètres, ou composantes du modèle.

Pour mettre en œuvre 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 points dans la grille de sortie (énergies £")■ ^{Pou r} i = l, - , n et ; = 1, ... , ^η', on dispose des données caractéristiques des intensités spectrales calculées À_i;- pour une énergie d'entrée E_t (données définies par l'utilisateur) et une énergie de sortie E (données mesurées). Pour illustrer le procédé selon l'invention, on suppose que les valeurs À_i;- ont été obtenues par l'intermédiaire d'un logiciel de simulation Monte-Carlo. Par conséquent, les données des intensités spectrales À_i;- sont considérées comme des réalisations issues d'une distribution de Poisson dont l'intensité sera estimée au moyen d'une procédure de régression non paramétrique. Comme Xi_j est alors une quantité positive ou nulle, on introduit y_tj = log(A_i;- + ε) où 0 < ε « 1. Dans ces conditions, la distribution de probabilité de y_tj peut être approximée pour des valeurs suffisamment grandes de (typiquement >

10) par une loi gaussienne de moyenne log l_£,- et de variance .

La distribution de poisson sera exprimée de la manière suivante : Distribution Poisson(l), pour n G N,

_ λ^η __{λ η}

/Poisson (λ) ( ) -—_\ ^{e n}

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, E y) à partir des points simulés (Ει, Ε ,γ^), afin d'en déduire le modèle 5(£^", £") pour tout (£, £") E KL² où IL est un espace continu :

où le symbole E(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 f{E, E , y), on choisit, par exemple, un mélange par processus de Dirichlet (DPM) comme distribution a priori. Cette modélisation s'exprime à partir d'une distribution aléatoire G telle que

G ~ DPO, G₀)

où le symbole "~" signifie "est distribué suivant" et DP(a, G₀) 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, Dunson D B and Quintana F 2010 Bayesian Nonparametrics Cambridge Séries in Statistical and Probabilistic Mathematics) avec une mesure moyenne G₀ 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 du DP est issue de Sethuraman (1994) qui exprime la mesure aléatoire G(-) comme :

avec w₁ = V_lt w_k = V_k a) et 5_U(-) représente la fonction de Dirac localisée en u. Ici, 9_k~G₀ représentent les paramètres associés à la k^ième composante de G. Le caractère non paramétrique est issu du nombre de composante potentiellement infini intervenant dans la somme caractérisant l 'a priori sur G.

La distribution D\ùch\et( ₁, ₂, ... , _K) est définie avec x = (x₁,x₂, ... , x_K)^' e R^K, xi > 0 pour 1 < i < K, E -To¹ ^ < 1, x_K = 1 -∑ϊ^₀ ^{1 χ}ί ^et i > ⁰ P^{ou r} 1≤ ί≤ K, f r_v-\ _ Γ(α₁+α₂+-·+α_κ) „, a₁-l_v a₂-l _v K-l

/Dirichlet(_ai,a₂ a_K) W - _Γ(_αι)Γ(α₂)...Γ(α ) ^{Xl 2} -^Χκ '

La distribution Beta(a, b), est définie pour 0 < x < 1 :

La distribution aléatoire f(E,E,y) s'exprimera alors : f(E,E^',y) = ^w_kfo_k{E,E^',y)

k=l

où fe_k indique la k^ième composante de f{E,E,y), paramétrée par 9_k.

Les composantes f_Qk 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 :

• 0_k = G½ A, a ec _¾ = (μ¾,μ'_¾), r_k = (τ_¾,τ'_¾),

• X_K(É) le vecteur de régresseurs centré avec Ë = (Ε,Ε^'), et _k le vecteur de coefficients de régression. On remarquera que le modèle de régression pour l'issue y sachant les covariables (Ε,Ε^') peut être choisi linéaire (i.e.

T

X_k{É) = (l,E - μ^Ε^' - μ'_¾) où v^T représente la transposée du vecteur v) ou polynomial dans l'approche DPGLM (ex. X_K(Ë) = (l,E - μ^Ε - -μ_¾)(£ -μ'_¾),(£ - μ_κ ² , (E - μ' _k) ) pour une régression quadratique),

• De plus, on définit la matrice ∑_fe, dépendant du paramètre xfj_k G {0,1}, comme suit: ∑_k = R^.^ _Τ°)-^ΒΨ* ^{AVEC R} _k = (J °) ^si ^ = ^{0 et}

^RPk ⁼ ~r(i ^~i ) ^si ^fe ^{= La variaDle} discrète binaire permet ainsi de choisir entre une composante alignée sur les axes E et E {ifj_k = 0) et une composante oblique orientée par la droite E = E (4>_k = 1).

Munis de ces notations, chaque composante de la distribution aléatoire s'exprime comme : fe_k(E, E^',y) = N₂

où, Χ · \μ, σ²) représente la loi gaussienne de μ et de variance σ² , Κ₂ {· \μ, Τ) la loi gaussienne bivariée de moyenne β G E² et de matrice de covariance∑.

Une distribution gaussienne multivariée Ν_ρ(μ, Ω) est définie pour X G R^p,

1 î

où |4| représente ici le déterminant de la matrice Note : cette définition couvre également le cas univarié quand p = 1.

La mesure a pr/or/ G₀ du processus de Dirichlet est choisie de la façon suivante. La loi a priori du paramètre _¾ est une gaussienne, la variance r_k est distribuée suivant une loi gamma-inverse, ifj_k suit une loi a priori de type Bernoulli et l'a priori pour les coefficients de régression coefficients β_ίι est une loi normale (gaussienne) multivariée de dimension

En utilisant la règle de Bayes, nous déduisons de l'expression de f(E, E^',y) que : f(y ï) =∑ _{y Λ} » (r|ol^■ *■(¾)■«- *'">)■ Finalement, en calculant l'espérance conditionnelle E E^') à partir de l'expression de E^'), on obtient l'expression du modèle de réponse s ectrale S(E, E^') :

Cette dernière expression permet de rendre compte de l'a priori correspondant à l'approche bayésienne non paramétrique pour la réponse spectrale débruitée et interpolée. Les caractères hétéroscédastiques et non gaussiens en tout (Ε, Ε^') sont pris en compte dans le modèle de E^') grâce au mélange de distributions gaussiennes M ( |/?[ · X_K (É), _E

A partir de ce choix de modèle probabiliste et des données observées par simulation physique (E_£ le procédé va estimer la loi a posteriori f{E, ... , Ε_η, Ε'_η;γ_ηη·) et l'espérance conditionnelle de cette loi a posteriori : S(E, E^') = E yE, E , E_x, E , y_n, ... , E_n, E'_n; y_nn< représentant l'intensité dans le spectre d'une énergie de sortie E , pour toute énergie d'entrée E.

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é adopte une approche dite d'échantillonnage par tranche d'après Kalli et al. (Kalli M, Griffin J E and Walker SG 201 1 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 Pour T générations, la distribution a posteriori de la réponse spectrale est approximée par l'ensemble des tirages S(E, E ^')^ pour t = l, ... , T, et la réponse estimée (moyenne a posteriori) s'exprime :

Il est également possible de calculer de façon analogue l'écart-type a posteriori et les intervalles crédibles à partir de la collection

Un exemple pour exécuter la procédure d'échantillonnage selon Gibbs va maintenant être donné. Le modèle génératif pour les paramètres du mélange DPM est donné par :

V_X,V ,... ~ Beta(l,a)

soient w₁ = V_ltw₂ = y₂(l - Vi),■■■

τι ¹> ^τ2¹>■■■ ~ Gamma(a_T, b_T)

τ^,τ^¹, Gamma(a_T.,ô_T.)

β₁,β₂,...~Μ_ι^₎₁{Μ_β,Γ_β)

-φ₁,-φ₂> -~ Bernoulli(p)

soient θ₁ = (β^τ^β^-ψ^, θ₂ = (μ₂,τ₂,β₂,ψ₂), ...

Ce modèle génératif est entièrement équivalent à la densité de probabilité a priori f(w₁,w₂, ...,θ_χ, θ₂, ...), et correspond à un a priori de type processus de Dirichlet pour la mesure aléatoire G(-) =∑_¾=1w_fe ¾(·) ~ ÛP a,G_Q).

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 Pour cela, on introduit des variables latentes dites de classification K^, définies pour i = l,...,n et j = l,...,n', tel que Ki = k si (Ei.E'j.yij) est distribué suivant la k^ième composante du mélange 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≤ n,

avec μ_Κί. = (μ_Κί μ'_κ.),∑*_y = , τ'*_.. "^ y_ij\ Κ_φΕ₀Ε), θ,, θ₂, ... ~ M (y_y I /¾^■ X_Ki]{E_itE ) ^^K^ ⁱ ). Ce modèle génératif est entièrement équivalent à la vraisemblance :

f(.Kn_>■■■ , Knn'> E±, E i, ... , E_n, E _n;

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 θ₂, ... I K ₁, ... ,K_nn',E_x,E _x, ...,E_n,E _n',yii, — ,y_nn'- 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 suivants, pour tout k,

w_1}w₂, ... | K_11} ...,Κ_ηη·

I fej k> Kn>^■■■ > Knn'> Ei, E -L, ... , E_n, E _η·

^Tk \ E _x, ... , E_n, E _η· τ k \ fej k'^k' Kll> — , Κηη'Έί, Ε ... ,E_n,E _n ^> fej fc E _x, ... , E_n, E _n ^> k \ fej k'^Tk' ... , E_n, E _n',yii, — ,y_nn' 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≤ n,

ΐ μ_χ,μ₂, ... ,μ _χ,μ _v ... ,τ_χ,τ₂, ... ,τ _χ,τ ₂, ... ,β_χ,β₂, ...,ψ₁,ψ₂, ... ,w₁,w₂, ...

Le comportement non paramétrique du procédé est caractérisé par la collection potentiellement infinie de paramètres w_k et 9_k. Afin de 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 aléatoire fini κ 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 : Initialisation aléatoire.

- Générer pour 1 < i≤ n,j≤ n',u_i;~Uniform(0, 1/n.n'), et poser u^* = min_u({ y})

- Générer l^-BetaCl.a) ; Poser w₁ = V₁ and r_x = 1 - V.

- Pour k > 1, générer V_k~Beta(l,c), poser w_k =V_kr_k_₁, poser r_k =

r^U-l^tant que r_k≥u^*.

- Affecter à K* la valeur maximale de k.

- Générer 6_k pour 1 < k≤ K*, à partir des distributions a priori.

À l'itération t = l,...,T,

· Générer

w₂, ...) pour ί≤ η, j≤ η'

- Calculer pour 1 < k≤ κ^*,

- X_k = (w_k > u_y)max(w_k^/n.nD:Ar₂(£_i;|^

1(A) est la fonction indicatrice : 1(A) = 1 si A est vrai et 1(A) = 0

sinon.

- Générer Κ_ίΓ -J—∑ i

Lk=\^Ak

• Réordonner les étiquettes de composantes en suivant leur ordre

d'apparition dans la génération des Κ^. Affecter à κ_η le nombre de valeurs distinctes de Κ^. Poser, pour tout k≤ κ_η, n_k = #{κ^ = k], le nombre d'observations y_i;- assigné à la composante k.

• Générer (w_1; w₂... , u , ... , u_nn\ K_llt ... , K_nn

- Générer (w_lt w₂, ... , w_Kn, r_Kn) ~ Dirichlet(n₁, n₂, ... , n_Kn, a)

- Générer pour i < n,j≤ n', it_i;~Uniform(0,min(w_/c, 1/n.n'))

- Poser u^* = min_u({tt_i;})

- Pour k > κ_η, générer l^-BetaC^a), poser w_k = V_kr_k_₁, poser r_k = r_fe--1(l - V_k), tant que r_k≥ u^*.

- Affecter à K* la valeur maximale de k.

• Générer ...,E_n,E'_n) pour k≤ K*,

Calculer

- Générer ν ~ Uniform(0,1)

- Si <—— , poser ip_k = 0 sinon poser x/j_k = 1

Vo+Vi

• Générer (τ_¾| μ^μ'^ψ^Κ^, K_nn >,E_XlE ,...,E_n,E'_n) pour /c < /^*,

- Si = 0 T_fc ¹~Gamma

- = l -μ',)

· Générer (τ^\μ^μ'_]ί,ψ_1{,Κ₁₁,...,Κ_ηη;Ε₁,Ε'₁,...,Ε_η,Ε'_η·) pour/c < κ^*

- Si = 0

-

• Générer (ji_k \ T_k,r'_k, ip_k, K_llt ... , K_nn; E_x, E , ... , E_n, E'_n) pour k≤ κ

- Calculer = (∑_{ij. _K._j=k] E_t ,∑_{[ij: K}._j=k} E ) - Poser ∑_k = R^_k.^ _T°)- _k avec _k = (J J) si tf_fc = 0 et - Calculer ν_μ = (Γ_μ ^_1 + n_k∑_fe ^_1)^~

- Calculer Μ_μ¾ = ν_μκ. (Γ_μ ^_1.τ^η _μ +∑_fc ^_1-¾

- Générer μ^- ₂{Μ_μκ,ν_μκ)

• Générer μ^μ'_ίι>^τ^^τ >'Φ^^ιι_> ... ,Κ_ηη<,Ε_χ,Ε , ...,Ε_η,Ε'_η<,γ_χχ, ...,y_nn) pour k≤ K*,

- Poser

- Poser

- Générer ^-^^(Μ^,Γ^).

Pour ce tirage aléatoire, on approxime la variance d'observation de yt_j (elle-même approximée à une variable aléatoire gaussienne) σ?- = _e-fiI-Xk(Èⁱj) à la valeur de l'observation Poisson A_i;- = e^~yi) correspondante. Cette approximation, qui simplifie grandement l'inférence, est d'autant plus valable que A_£y- est élevé.

• Calculer à chaque itération (t) pour tout choix de grille de prédiction Ë = (E,E'),

La réponse spectrale débruitée peut, au choix, partager ou non la même grille (E_£,E ) que la simulation Monte-Carlo physique initiale. En effet, l'utilisateur peut choisir n'importe quel point (£,£") d'intérêt. En particulier, il est possible d'interpoler entre les points de la grille initiale.

Finalement, à partir de T itérations de l'échantillonneur de Gibbs, on obtient une estimation du méta-modèle pour la spectrométrie τ

S ÇE. E') * - } S ÇE. E^

t=i

Il est alors possible de calculer l'écart-type a posteriori de l'estimateur de réponse spectrale

t=i

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 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(<p_fc, ¾) sur les paramètres d'échelle b_T et b dans la distribution des amplitudes de composantes. Ce choix conduit à une distribution a posteriori Gamma.

Exemple :

• 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).

« Le paramètre p de la loi de Bernoulli intervenant dans le modèle de mélange du type de composante ifj_k peut être lui-même estimé par l'intermédiaire d'un hyperprior BetaC^, ^)- La loi a posteriori pour p suit alors l'expression :

p~ Beta + #{k≤ κ_η ·. i _k = 0}, _πι + #{k≤ κ_η ·. i _k = 1})

• Notons également qu'il peut être efficace de remplacer Va priori gaussien sur les coefficients de régression _k par un a priori dit empirique basé sur le tirage aléatoire de P points { ^ Έ'^ ... , Ë_P, Ë'_P) à partir de ·7\Γ₂( _¾, Σ_¾) et de prendre pour ( _1; ... , y_P) la plus proche valeur de y_tj correspondant à chaque point échantillonné. P doit être plus grand que la taille du vecteur k. Nous générons alors β_1ί ~ Ν(Μ_β, Γ_β) comme loi a priori (avec

Ep (^p E ) )■

L'expression pour la loi a posteriori de _k 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 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 nuclear spectrometry," Nuclear Science Symposium Conférence Record, 2007. NSS '07. I EEE , 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 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 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. 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 mono-énergétiques, un a priori par processus de Dirichlet sur les radioéléments d'une base de données de radionucléïdes. 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 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 radionucléïdes 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 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éïdes, on utilise un méta-modèle étendu noté Ξ(Ε, Ε^', ξ) où ξ représente un paramètre caractéristique de l'effet de matrice associé à la source. La collection d'effets 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.

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 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 χ 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 élément par un indice entier. On note N_x le nombre d'émetteurs retenus pour l'élément χ considéré et n_x pour l = 1, ... , N_X, les probabilités d'émissions associées ainsi que v_xX pour l = 1, ... , N_X, les énergies correspondantes. On définit la réponse du radionucléide ψ_χ(Ε', ξ), pour une énergie observée £" et un effet de matrice ξ, par :

on définit Α_{χ ξ} = /_ο ^°° t l = 1, ... , N_X. On définit

alors, une réponse de radionucléide normalisée par :

1 = 1

on vérifie que J₀ ^°° i/ ^* _ξ (£") dE' = 1.

La densité de probabilité du i^eme photon observé dans le canal d'énergie s'exprime alors pour i = 1, ... , n:

k= l

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 μ_ρ et d'écart-type σ_ρ :

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, F°XF^) Ρχ{χ) ^sj(x) représente la distribution de probabilité a priori d'observer un élément χ où / est le nombre de radionucléides de la base retenue pour l'analyse et p° > 0. 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 du processus :

fc=l

avec w₁ = V_x, w_k = et 5_U(-) 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 K_t d'allocation du i^eme photon 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,

...,χ_ί,χ₂, ...,ξ_ί,ξ₂, ...,ρ

et pour tout k

w₁,w₂, ... | K_lt ...,K_n

X_k \ Κ-_L, ...,K_n,E_lt ... , Ε_η,ξ₁,ξ₂, ...,p

Κι_>—,K_n,E , -,Ε_η,χ₁,χ₂, ...,p

ainsi que

p| K_lt ...,K_n,E_lt■■·,Ε_η,χ₁,χ₂, ...,ξι,ξ₂,■■■

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 _£ qui permet de ne générer qu'un nombre aléatoire fini κ de composantes à l'itération (t) tout en évitant une troncature arbitraire du modèle sont utilisées. Sous ces hypothèses, les différentes étapes de génération aléatoire de l'algorithme sont détaillées ci- après.

Initialisation aléatoire.

- Générer pour 1 < i < n, ιι_έ~11ηιτοπτι(0, 1/n), et poser u^* = min_u({uj) - Générer l^-BetaCl, a) ; poser w₁ = V₁ et r_x = l- V .

- Pour k > 1, générer V_k~Beta(l, a) , poser w_fe = V_fer_fe_₁, et poser r_k = r_fc__!(l - Vfc), tant que r_fe > u^*.

- Affecter à K* la valeur maximale de k.

- Générer χ_κ pour 1 < /c < K*, à partir de la distribution a priori F°.

- Générer ξ_}ι pour 1 < k≤ K*, à partir de la distribution a priori F°.

- Générer p ~

À l'itération t = l,...,T,

• Générer (K_I\E_i ^',w₁,w₂, ...,x_ltx₂, ...,ξ₁,ξ₂, ...,p) pouri < n,

- Calculer pour 1 < k≤ K*,

- À_k = l(w_k >u_t) max(w_fe,l/n) ^^* _¾¾( £^"-) où 1(A) est la fonction indicatrice : 1(A) = 1 si A est vrai et 1(A) = 0 sinon.

- Générer K_T~ -J—∑ i S_K(KÔ,

• Réordonner les étiquettes de composantes en suivant leur ordre d'apparition dans la génération des K_t. Affecter à κ_η le nombre de valeurs distinctes de K_t. Poser, pour tout k < κ_η, n_k = #{K_t = k}, le nombre d'observations assignées à la composante k.

• Générer (w₁,w₂ ...,u_n\ K , ...,K_n)

- Générer (w_lt w₂, ... , w_Kn, r_KJ ~ Dirichlet(rr_lJ n₂, ... , n_Kn, a),

- Générer pour i < 1/n)),

- Poser u^* = min_u({uj),

- Pour k > K_n, générer l^-BetaC^a), poser w_k = V_kr_k_₁, poser r_k = r_fe--1(l - V_k), tant que r_k≥ u^*,

- Affecter à K* la valeur maximale de k.

· Générer Κ₁,...,Κ_η,Ε₁ ^',...,Ε_η,ξ₁,ξ₂,...,ρ) pour k≤ K*,

- Calculer pour tout; = 1, ...,],

- Générer^- —∑^J _j=1Vk (x_k).

Lj₌₁Vk,j Générer (ξ \ Κ_1} ... , Κ_η, Ε , ... , Ε_η,χ₁,χ₂, ... , p) pour k

- Calculer pour tout m =

• Générer p\ ...,K_n,E , ...,Ε_η,χ₁,χ₂, ...,,ξ₁,ξ₂, ...) par une étape Metropolis-Hastings,

- Générer une marche aléatoire d'écart-type ε_ρ autour de p

- Calculer

- Générer une variable uniforme v ~ Uniform(0,1).

- Si v≤ min l,g_MH) 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,

De surcroît, il est également immédiat de calculer l'écart-type a posteriori des activités :

Il est également possible d'obtenir le spectre d'entrée estimé, déconvolué de la réponse du système

N (t)

T %k

^® ^¾ r ^wit) ¾^ ¾^(f)

t=i 1=1 ^K 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 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 l'invention. La courbe en pointillé représente le spectre simulé, la courbe en trait plein le spectre reconstruit.

Les figures 1 1 et 12 illustrent un deuxième exemple de spectre. La figure 1 1 représente le spectre des énergies incidentes, avec une région d'émission du premier pic au ⁶⁰Co, 1 173,2 keV. La figure 12 le même spectre avec un zoom sur la région d'intérêt.

Claims

REVENDICATIONS

1 - 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,

· 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 £", • pour i = l, ... , n et j = l, ... , n', on dispose des données caractéristiques des intensités spectrales calculées À_£y- pour une énergie d'entrée £_£ et une énergie de sortie £",-,

• 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 (Ει, Ε , γ^), et on déduit un modèle S(E, E') pour tout (E, E') E R² OÙ M est un espace continu :

où le symbole E(y) représente l'espérance mathématique de la variable aléatoire y, y est déduite des intensités spectrales calculées λ ·, • à 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 - Procédé selon la revendication 1 caractérisé en ce que l'on détermine les valeurs À_£y- en utilisant un logiciel Monte-Carlo, les données À_£y- étant considérées comme des réalisations issues d'une distribution de Poisson dont l'intensité est estimée au moyen d'une procédure de régression non paramétrique et on introduit y_tj = log(A_i;- + ε) où 0 < ε « 1, puis

on approxime la distribution de probabilité de y_tj pour des valeurs suffisamment grandes de > 10 par une loi gaussienne de moyenne log et de variance , pour la loi jointe f(E, E' ,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 f_Qk composantes de f(E, E ,y), f{E, E^', y) = ^ w_k fo_k{E, E^',y)

k=l

paramétrée par 9_k le paramètre associé à la k^ième composante de G mesure aléatoire définie par

G(-) =∑w_fc ¾(·)

fc=l

avec w₁ = V_x, w_k = Ι^ Π?=ι (1 - ½) tel que V_k~Be\.a(l, a) et 5_U(-) représente la fonction de Dirac localisée en u et Beta(a, ô), pour 0 < x < 1 avec

3 - Procédé selon la revendication 2 caractérisé en ce que les composantes de la loi jointe f_Qk sont exprimées à partir des valeurs suivantes :

• 0_k = G½ A, a ec _¾ = (μ¾, μ'_¾), r_k = (τ_¾, τ'_¾) , • X_K(Ê) le vecteur de régresseurs centré avec Ê = (Ε,Ε^'), et _k le vecteur de coefficients de régression,

• la matrice ∑_fe, dépendant du paramètre ifj_k G {0,1}, comme suit : ∑_k = xfj_k = 1, permettant de choisir entre une composante alignée sur les axes E et E (4>_k = 0) et une composante oblique orientée par la droite E = E W>_k = !)^■

chaque composante _θ¾ s'exprime alors :

où, ^\ (·|μ,σ²) représente la loi gaussienne de μ et de variance σ², Κ₂{· \μ,Τ) la loi gaussienne bivariée de moyenne μ G E² et de matrice de covariance∑

où

la loi a priori du paramètre _¾ est une gaussienne, la variance r_k est distribuée suivant une loi gamma-inverse, ifj_k suit une loi a priori de type Bernoulli et l'a priori pour les coefficients de régression coefficients _k est une loi normale (gaussienne) multivariée de dimension et

en appliquant une règle de Bayes sur l'expression de f{E,E,y) on obtient l'expression

et le modèle probabiliste A partir de ce modèle probabiliste et des données observées par simulation {EuE'_j.yi_j) , on estime la loi a posteriori f E ,E',y\E₁,E'₁,y₁₁, ...,E_n,E'_n y_nni) et l'espérance conditionnelle S E,E') = E y\E,E' ' ,E_X,E' ,y_xl, ...,E_n,E'_n,,y_nni) afin de déterminer les éléments présents dans l'objet et leur activité. 4 - Procédé selon la revendication 3 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ératio procédure MCMC, on génère une réponse spectrale débruitée

pour T générations, la distribution a posteriori de la réponse spectrale est approximée par l'ensemble des tirages S^E. E ") ^ pour t = l, ... , T, et la réponse estimée s'exprime :

T

S{E, E^') ~ - S(E, E^')

t=i

5 - Procédé selon la revendication 4 caractérisé en ce que le schéma d'approximation comporte une étape d'échantillonnage par tranche utilisant un nombre aléatoire fini κ de composantes pour chaque itération et en ce qu'il comporte les étapes suivantes :

on introduit des variables latentes de classification K^, définies pour ί =

Ι, .,. , η et ; = Ι, .,. , η', tel que K_t = k si est distribué suivant la fciè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 < n ,

équivalent à la vraisemblance (#11,■■■ , Knn'_> E±, E i, ... , E_n, E _n; yn_> -_>y_nn'\ w₁,w₂, ...,Θ_χ, e₂, ...) à partir de ces distributions de probabilités et par application de la règle de Bayes, on calcule la densité de probabilité conditionnelle

f(w₁,w₂, ...,θ , θ₂, ... I fin, ...,K_nn; E_x, E , ...,E_n,E'_n>,y_xx, ...,y_nn

en utilisant un échantillonner de Gibbs qui à chaque itération (t) génère successivement les échantillons suivants, pour tout k,

w_1}w₂, ... I #11, ...,#_ηη·

I fej k> #11 > Knn'> Ei> E i,^■■■ , E_n, E _n ^>

Tk \ _x, ... , E_n, E _η· τ k \ _x, ... ,E_n,E _n ^> fej E _x, ... , E_n, E _n ^> k \ fej k'^Tk ... , E_n, E _n', n, ... ,y_nn' 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≤ n,

ΐ μ_χ,μ₂, ... ,μ _χ,μ _v ... ,τ_χ,τ₂, ... ,τ _χ,τ ₂, ... ,β_χ,β₂, ...,ψ₁,ψ₂, ... ,w₁,w₂, ... à partir de T itérations de l'échantillonneur de Gibbs, on obtient une estimation du méta-modèle pour la spectrométrie : 6 - Procédé selon la revendication 5 caractérisé en ce que l'on introduit une variable auxiliaire afin de ne générer qu'un nombre aléatoire fini κ de composantes à l'itération (t) tout en évitant une troncature arbitraire du modèle. 7 - Procédé selon l'une des revendications 4 à 6 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')^}

8 - Procédé selon l'une des revendications 1 à 7 caractérisé en ce que l'on introduit un a priori Gamma(<p_fc, ¾) sur les paramètres d'échelle b_T et b'_T dans la distribution des amplitudes de composantes conduisant à une distribution a posteriori Gamma :

avec Gamma(a, ô), pour > 0,

/Gamma(a,6) (·*·) ^— Γ(α) ^

9 - Procédé selon la revendication 3 caractérisé en ce que l'on remplace l 'a priori gaussien sur les coefficients de régression _k par un a priori basé sur le tirage aléatoire de P points { ^ Ε'^ ... , Ë_P, È'_P) à partir de J\f₂ Q½,∑_fe) et on prend pour ( _1; ... , y_P) la plus proche valeur de y_tj correspondant à chaque point échantillonné où P est plus grand que la taille du vecteur _k, on génère alors _k ~ Ν(Μ_β, Γ_β) comme loi a priori (avec Ë = (Ë_p, Ë'_p))

10 - 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é Ξ(Ε, Ε^', ξ) où ξ est un paramètre identifié par un indice entier et caractéristique d'un effet de matrice. 1 1 - 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 Ν_χ le nombre d'émetteurs retenus pour l'élément χ considéré et π_χ pour 1 = 1, ... , N_X, les probabilités d'émissions associées ainsi que v_x pour l = 1, ... , Ν_χ, les énergies correspondantes,

on définit la réponse du radionucléide -ψ_χ{Ε', ξ), pour une énergie observée E' et un effet de matrice ξ, par

on définit .^ = J₀ ^°° ^^(£") d£" et π* _{ι ξ} = pour tout / = 1, ... , N_X, on définit une réponse de radionucléide normalisée par

1=1

on vérifie que J₀ ¾(£") dE' = 1

on exprime la densité de probabilité du i^eme photon observé dans le canal d'énergie pour i = Ι, .,. , η f{E[) = _jWk r_Xk,_?k{p- ^Ei')

k=l

où p est un paramètre positif de conversion de l'indice de canal en énergie (keV/canal).

On introduit les variables K_t d'allocation du i^eme photon observé à une composante k du mélange,

On alterne des tirages aléatoires suivant des lois conditionnelles suivantes, pour tout i≤ n,

pour tout k

w₁, w₂, ... | K_lt ... , K_n Xk \ K-L, ...,K_n,E_lt ... , Ε_η,ξ₁,ξ₂, ...,p

Κι_>—,K_n,E , -,Ε_η,χ₁,χ₂, ...,p

ainsi que

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

5 k,

t=l

avec p un a priori gaussien centré sur μ_ρ et d'écart-type σ_ρ.

12 - Procédé selon la revendication 11 caractérisé en ce que l'on calcule o l'écart-type a posteriori des activités

et/ou le spectre d'entrée estimé, déconvolué de la ré onse du système

5