Étude Approfondie : Détection Automatique des Isotopes Radioactifs avec un Modèle d'Apprentissage Machine en Python

Analyse Exploratoire des Données : Spectroscopie Gamma en Python (Partie 2) Dans cette deuxième partie, nous allons approfondir l'analyse exploratoire des données de spectroscopie gamma en montrant comment créer et entraîner un modèle d'apprentissage automatique (machine learning) pour détecter les éléments radioactifs. Préalable Important Avant de commencer, un avertissement crucial : toutes les fichiers de données collectés pour cet article sont disponibles sur Kaggle, ce qui permet aux lecteurs de former et de tester leurs modèles ML sans disposer d'un matériel réel. Si vous souhaitez tester des objets réels, faites-le à vos risques et périls. J'ai utilisé des sources d'échantillons qui peuvent être légalement trouvées et achetées, comme du verre d'uranium antique ou des montres anciennes à peinture de radium. Vérifiez toujours les lois locales et suivez les directives de sécurité relatives à la manipulation de matériaux radioactifs. Bien que ces sources ne soient pas très dangereuses, elles doivent être manipulées avec précaution. Objectif L'objectif est de créer un modèle ML capable de détecter automatiquement divers éléments radioactifs en utilisant un détecteur de scintillation Radiacode. Méthodologie 1. Spectre Gamma Un spectre gamma est une trace des énergies des rayons gamma émis par un objet radioactif. Contrairement au compteur Geiger qui ne mesure que le nombre de particules radioactives, un détecteur de scintillation indique également les énergies de ces particules. Cela fournit des informations essentielles sur les isotopes présents. Pendentif Ionogène : Acheté dans une boutique chinoise, ce pendentif annoncé comme « générateur d'ions » présente un niveau de radioactivité de 1,20 µSv/h, soit 12 fois celui du fond normal (0,1 µSv/h). Son spectre gamma révèle des isotopes de thoriure-232. Minéral (Uraninite) : Ce spécimen peut être trouvé dans certaines régions comme l'Allemagne, la République Tchèque ou les États-Unis. Son spectre gamma montre des pics caractéristiques de l'uranium mais pas de thoriure, ce qui permet de distinguer les éléments. Collecte des Données La collecte des échantillons et des données constitue un défi majeur, notamment en raison des restrictions légales et des mesures de sécurité nécessaires. Les échantillons choisis incluent des détecteurs de fumée à l'américium, des montres anciennes à peinture de radium, du verre d'uranium et des baguettes de tungstène thoriated. La collecte de chaque spectre prend 10 à 20 minutes, et 10 à 20 échantillons par objet sont recommandés, ce qui peut demander plusieurs heures. Au total, 100 à 200 fichiers JSON peuvent être collectés. Le script Python utilisé pour collecter et sauvegarder les spectres est le suivant : ```python from radiacode import RadiaCode, RawData, Spectrum import random import time import json from datetime import datetime import logging def read_forever(rc: RadiaCode): """ Lecture continuelle des données du détecteur """ while True: interval_sec = random.randint(1060, 3060) read_spectrum(rc, interval_sec) def read_spectrum(rc: RadiaCode, interval: int): """ Collecte et sauvegarde du spectre """ rc.spectrum_reset() dt = datetime.now() filename = dt.strftime("spectrum-%Y%m%d%H%M%S.json") logging.debug(f"Collecte du spectre pour {interval // 60} min") t_start = time.monotonic() while time.monotonic() - t_start < interval: show_device_data(rc) time.sleep(0.4) spectrum = rc.spectrum() spectrum_save(spectrum, filename) def show_device_data(rc: RadiaCode): """ Affiche les taux de décompte (CPS) """ data = rc.data_buf() for record in data: if isinstance(record, RawData): logging.debug(f"CPS: {int(record.count_rate)}") def spectrum_save(spectrum: Spectrum, filename: str): """ Sauvegarde les données du spectre """ duration_sec = spectrum.duration.total_seconds() data = { "a0": spectrum.a0, "a1": spectrum.a1, "a2": spectrum.a2, "counts": spectrum.counts, "duration": duration_sec, } with open(filename, "w") as f_out: json.dump(data, f_out, indent=4) logging.debug(f"Fichier '{filename}' sauvegardé") rc = RadiaCode() read_forever(rc) ``` Préparation des Données pour l'Entraînement du Modèle 2. Chargement des Données Les spectres gamma sont stockés au format JSON. Un fichier JSON individuel ressemble à ceci : json { "a0": 24.524023056030273, "a1": 2.2699732780456543, "a2": 0.0004327862989157, "counts": [48, 52, ..., 0, 35], "duration": 1364.0 } Nous chargeons ces fichiers avec Python : ```python import json from dataclasses import dataclass @dataclass class Spectrum: duration: int a0: float a1: float a2: float counts: list[int] def channel_to_energy(self, ch: int) -> float: return self.a0 + self.a1 * ch + self.a2 * ch**2 def load_spectrum_json(filename: str) -> Spectrum: with open(filename) as f_in: data = json.load(f_in) return Spectrum( a0=data["a0"], a1=data["a1"], a2=data["a2"], counts=data["counts"], duration=int(data["duration"]), ) ``` 3. Transformation des Données Pour uniformiser les données, je les normalise et ajoute du bruit pour augmenter la diversité. 3.1. Normalisation Le nombre de particules détectées peut varier selon la durée de collecte et l'intensité de la source. Pour résoudre ce problème, j'utilise une filtre de Savitzky-Golay pour réduire le bruit et normalise ensuite les données entre 0 et 1. ```python import numpy as np from scipy.signal import savgol_filter def smooth_data(data: np.array) -> np.array: window_size = 10 data_out = savgol_filter( data, window_length=window_size, polyorder=2, ) return np.clip(data_out, a_min=0, a_max=None) def normalize(spectrum: Spectrum) -> Spectrum: counts = np.array(spectrum.counts).astype(np.float64) counts = smooth_data(counts) val_norm = counts.max() return Spectrum( duration=spectrum.duration, a0=spectrum.a0, a1=spectrum.a1, a2=spectrum.a2, counts=counts / val_norm ) ``` 3.2. Augmentation des Données Pour augmenter le nombre de spectres, j'ajoute du bruit aux échantillons originaux. Je sélectionne le canal 680 keV, où il n'y a pas d'isotopes intéressants, et ajoute du bruit à 50% de l'amplitude de ce canal. python def add_noise(spectrum: Spectrum) -> Spectrum: counts = np.array(spectrum.counts) ch_empty = spectrum.energy_to_channel(680.0) val_norm = counts[ch_empty] ampl = val_norm / 2 noise = np.random.normal(0, ampl, counts.shape) data_out = np.clip(counts + noise, a_min=0, a_max=None) return Spectrum( duration=spectrum.duration, a0=spectrum.a0, a1=spectrum.a1, a2=spectrum.a2, counts=data_out.tolist() ) Entraînement du Modèle 3.3. Chargement et Préparation des Données Je combine les fichiers de données en un seul ensemble pour l'entraînement du modèle. ```python import glob from sklearn.preprocessing import LabelEncoder all_files = [ ("Americium", glob.glob("../data/train/americium.json")), ("Radium", glob.glob("../data/train/radium.json")), ("Thorium", glob.glob("../data/train/thorium.json")), ("Uranium Glass", glob.glob("../data/train/uraniumGlass.json")), ("Uranium Glaze", glob.glob("../data/train/uraniumGlaze.json")), ("Uraninite", glob.glob("../data/train/uraninite.json")), ("Background", glob.glob("../data/train/background*.json")), ] def prepare_data(augmentation: int) -> Tuple[np.array, np.array]: x, y = [], [] for name, files in all_files: for filename in files: sp = normalize(load_spectrum_json(filename)) for _ in range(augmentation): sp_out = add_noise(sp) x.append(get_features(sp_out)) y.append(name) return np.array(x), np.array(y) ``` 3.4. Extraction des Caractéristiques Seules les valeurs des énergies des isotopes seront utilisées comme entrées du modèle. ```python isotopes = [ ("Am-241", 59.5), ("K-40", 1460.0), ("Ra-226", 186.2), ("Pb-214", 242.0), ("Pb-214", 295.2), ("Pb-214", 351.9), ("Bi-214", 609.3), ("Bi-214", 1120.3), ("Bi-214", 1764.5), ("Pb-212", 238.6), ("Ac-228", 338.2), ("TI-208", 583.2), ("AC-228", 911.2), ("AC-228", 969.0), ("Th-234", 63.3), ("Th-231", 84.2), ("Th-234", 92.4), ("Th-234", 92.8), ("U-235", 143.8), ("U-235", 185.7), ("U-235", 205.3), ("Pa-234m", 766.4), ("Pa-234m", 1000.9), ] def get_features(spectrum: Spectrum) -> np.array: energies = [energy for _, energy in isotopes] data = [spectrum.counts[spectrum.energy_to_channel(energy)] for energy in energies] return np.array(data) ``` 3.5. Entraînement du Modèle J'utilise un modèle XGBoost avec GridSearchCV pour trouver les meilleurs paramètres. Le modèle est ensuite enregistré pour une utilisation ultérieure. ```python from xgboost import XGBClassifier from sklearn.model_selection import GridSearchCV bst = XGBClassifier(n_estimators=10, max_depth=2, learning_rate=1) clf = GridSearchCV( bst, { "max_depth": [1, 2, 3, 4], "n_estimators": range(2, 20), "learning_rate": [0.001, 0.01, 0.1, 1.0, 10.0] }, verbose=1, n_jobs=1, cv=3, ) clf.fit(X_train, y_train) print("Best score:", clf.best_score_) print("Best params:", clf.best_params_) Enregistrement du modèle le.classes_ = np.load("../models/V1/LabelEncoder.npy") isotopes_save("../models/V1/isotopes.json") bst.save_model("../models/V1/XGBClassifier.json") np.save("../models/V1/LabelEncoder.npy", le.classes_) ``` Test du Modèle 4. Tests en Temps Réel Le code de test se connecte au détecteur Radiacode, lit les spectres gamma toutes les minutes et prédit les isotopes. ```python from radiacode import RadiaCode, RealTimeData, Spectrum import time import logging def read_spectrum(rc: RadiaCode): """ Lecture et prédiction du spectre """ spectrum = rc.spectrum() logging.debug(f"Spectre : {spectrum.duration} sec de collecte") result = predict_spectrum(spectrum) logging.debug(f"Prédiction : {result}") def predict_spectrum(spectrum: Spectrum) -> str: """ Prédiction des isotopes depuis le spectre """ features = get_features(normalize(spectrum)) preds = bst.predict([features]) return le.inverse_transform(preds)[0] def read_cps(rc: RadiaCode): """ Lecture des CPS (counts per second) """ data = rc.data_buf() for record in data: if isinstance(record, RealTimeData): logging.debug(f"CPS : {record.count_rate:.2f}") rc = RadiaCode() logging.basicConfig(level=logging.DEBUG, format="[%(asctime)-15s] %(message)s", datefmt="%Y-%m-%d %H:%M:%S") logging.debug("Modèle ML chargé") fw_version = rc.fw_version() logging.debug(f"Device connecté : firmware {fw_version[1]}") rc.spectrum_reset() while True: for _ in range(12): read_cps(rc) time.sleep(5.0) read_spectrum(rc) ``` J'ai testé le modèle avec un objet réel : une montre ancienne datant des années 1950 avec une peinture de radium. Son niveau de radioactivité est environ 5 fois plus élevé que le niveau foncier, mais reste sûr (environ 2 fois inférieur à celui d'un vol en avion). Conclusion Dans cet article, j'ai expliqué le processus de création d'un modèle d'apprentissage automatique pour prédire les isotopes radioactifs en utilisant un détecteur de scintillation Radiacode. J'ai aussi testé le modèle sur des échantillons radioactifs accessibles et légaux. Des voies d'amélioration possibles incluent : - Ajout de plus d'échantillons et d'isotopes. - Incorporation d'informations sur le niveau de radioactivité des objets. - Exploration d'autres types de modèles, comme les recherches vectorielles pour des résultats plus interprétables. Si vous êtes intéressé par ce sujet, le code complet est disponible sur ma page Patreon, et vous pouvez me suivre sur LinkedIn pour de plus petits posts périodiques. Merci de votre lecture. Profil de l'Entreprise Radiacode est un détecteur de radiation compact et abordable, parfait pour des expériences de spectroscopie gamma. Il est idéal pour les particuliers et les chercheurs indépendants souhaitant explorer la radioactivité sans avoir accès à des installations nucléaires de grande taille. Évaluation Professionnelle Selon des professionnels de l'industrie, la spectroscopie gamma demeure un domaine prometteur, en particulier avec l'avancée des technologies miniaturisées. Le modèle d'apprentissage automatique présenté ici offre une première étape solide vers une détection automatisée des isotopes radioactifs, bien que l'ajout de plus de données et de fonctionnalités pourrait le rendre encore plus précis et fiable.