Python으로 방사성 원소 감지 ML 모델 개발 및 테스트 방법 설명

감마 스펙트로스코피 데이터 분석: 파이썬을 활용한 탐색 (Part 2) 이 글에서는 감마 스펙트로스코피 데이터를 활용해 방사성 원소를 자동으로 감지하는 머신 러닝 모델을 만드는 방법을 설명하겠습니다. 이 프로젝트는 방사능 물질의 종류를 결정할 수 있는 중요한 단계입니다. 주의사항 데이터 수집에 사용된 모든 파일은 Kaggle에서 제공되며, 실제 하드웨어 없이도 모델을 훈련하고 테스트할 수 있습니다. 실제 물체를 테스트하고자 하는 독자는 반드시 현지 법률을 확인하고 방사능 물질 취급 안전 가이드라인을 읽어야 합니다. 본 테스트에 사용된 소스들은 법적으로 구매 가능한 것들로, 심각한 위험성이 없습니다. 그러나 여전히 주의해서 다뤄야 합니다. 감마 스펙트럼 이해하기 감마 스펙트럼은 방사능 물질이 발산하는 감마선의 에너지를 측정하는 것입니다. 이는 감마선이 다르게 발산되는 에너지 수준을 통해 물질의 종류를 알아낼 수 있게 해줍니다. 예를 들어, 이 중화국 상점에서 구입한 펜던트는 약간의 방사능을 가지고 있어 1,20 µSv/h의 방사능 수준을 보였습니다. 이는 배경 수준(0.1 µSv/h)의 12배지만, 비행기 탑승 시 받는 방사능 수준과 유사해 크게 위험하지 않습니다. 두 번째 예로, 우라니나이트라는 미네랄을 들 수 있습니다. 이는 많은 양의 산화우라늄을 함유하고 있으며, 독일, 체코, 또는 미국 일부 지역에서 발견될 수 있습니다. 감마 스펙트럼을 통해 이 미네랄이 우라늄을 함유하고 있지만 토륨은 포함하지 않는 것을 알 수 있습니다. 감마선은 실제로 포톤으로, 가시광선과 같은 전자기 스펙트럼에 속합니다. 따라서 방사능 물질은 각자의 고유한 '빛'을 내지만, 인간의 눈에는 보이지 않습니다. 또한, 최근 전자 기술의 발달로 감마 스펙트로스코피 장비를 중가형 스마트폰 가격으로 구입할 수 있게 되었습니다. 데이터 수집 첫 번째 도전 과제는 샘플을 수집하는 것입니다. 저는 핵 기관이 아니므로 채슘이나 스트론튬과 같은 교정된 시료에 접근할 수 없습니다. 하지만 아메리슘은 연기 감지기에, 라듐은 1960년대 이전 시계의 다이얼 페인트에, 우라늄은 1950년대 이전 유리 제조에, 토륨은 현재도 생산되는 토륨-텅스텐 막대를 쉽게 구할 수 있습니다. 자연 우라늄 광물도 미네랄 가게에서 구매 가능하지만, 안전 조치가 필요합니다. 감마 스펙트로스코피의 장점은 물체를 분해하거나 파손하지 않고 측정할 수 있다는 점입니다. 데이터 수집은 Radiacode 103G 감마선 검출기와 오픈소스 radiacode 라이브러리를 사용하여 수행되었습니다. 감마 스펙트럼 데이터는 Radiacode 안드로이드 앱을 통해 XML 형식으로 내보낼 수 있지만, 수작업은 시간이 많이 걸리고 번거롭습니다. 이를 해결하기 위해 Python 스크립트를 작성하여 임의의 시간 간격으로 데이터를 수집했습니다. ```python from radiacode import RadiaCode, RawData, Spectrum import random import time import datetime import logging import json from typing import List, Optional def read_forever(rc: RadiaCode): while True: interval_sec = random.randint(1060, 3060) read_spectrum(rc, interval_sec) def read_spectrum(rc: RadiaCode, interval: int): rc.spectrum_reset() dt = datetime.datetime.now() filename = dt.strftime("spectrum-%Y%m%d%H%M%S.json") logging.debug(f"Making spectrum for {interval // 60} min") t_start = time.monotonic() while time.monotonic() - t_start < interval: show_device_data(rc) time.sleep(0.4) spectrum: Spectrum = rc.spectrum() spectrum_save(spectrum, filename) def show_device_data(rc: RadiaCode): data = rc.data_buf() for record in data: if isinstance(record, RawData): log_str = f"CPS: {int(record.count_rate)}" logging.debug(log_str) def spectrum_save(spectrum: Spectrum, filename: str): 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"File '{filename}' saved") rc = RadiaCode() read_forever(rc) ``` 이 스크립트는 Radiacode 검출기를 물체 근처에 놓고 몇 시간 동안 실행하면 10-20개의 JSON 파일이 저장됩니다. 각 샘플에 대해 이 과정을 반복하면 최종적으로 100-200개의 파일을 수집할 수 있습니다. 모델 훈련 데이터 수집이 완료되면 모델 훈련을 시작할 수 있습니다. 모든 파일은 Kaggle에서 제공되므로 독자들도 자신의 모델을 만들 수 있습니다. 먼저, 데이터를 불러와 필요한 특징을 추출합니다. ```python import json import numpy as np from dataclasses import dataclass from typing import Optional @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 energy_to_channel(self, e: float): c = self.a0 - e return int( (np.sqrt(self.a1**2 - 4 * self.a2 * c) - self.a1) / (2 * self.a2) ) 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"]), ) ``` 이제 수집된 데이터를 시각화해 볼까요? ```python import matplotlib.pyplot as plt def draw_simple_spectrum(spectrum: Spectrum, title: Optional[str] = None): fig, ax = plt.subplots(figsize=(12, 3)) ax.spines["top"].set_color("lightgray") ax.spines["right"].set_color("lightgray") counts = spectrum.counts energy = [spectrum.channel_to_energy(x) for x in range(len(counts))] ax.bar(energy, counts, width=3.0, label="Counts") ticks_x = [spectrum.channel_to_energy(ch) for ch in range(0, len(counts), len(counts) // 20)] labels_x = [f"{ch:.1f}" for ch in ticks_x] ax.set_xticks(ticks_x, labels=labels_x) ax.set_xlim(energy[0], energy[-1]) plt.ylim(0, None) title_str = "감마 스펙트럼" if title is None else title ax.set_title(title_str) ax.set_xlabel("에너지, keV") plt.legend() fig.tight_layout() sp = load_spectrum_json("thorium-20250617012217.json") draw_simple_spectrum(sp) ``` 표준 지거 카운터는 방사능 입자의 수만을 알려주지만, 감마 스펙트럼을 이용하면 입자의 에너지 분포를 확인할 수 있습니다. 방사능 붕괴는 무작위적이므로, 측정 시간이 길수록 그래프가 부드러워집니다. 데이터 변환 노이즈 제거: 감마 스펙트럼은 짧은 측정 시간에 노이즈가 많이 발생할 수 있습니다. 이를 줄이기 위해 Savitzky-Golay 필터를 적용합니다. 정규화: 수집 시간과 소스 강도에 따라 검출된 입자의 수가 달라질 수 있으므로, 데이터를 0~1 범위로 정규화합니다. 데이터 증강: 샘플 수가 적으므로, 원래 스펙트럼에 노이즈를 추가해 더 많은 합성 샘플을 생성합니다. ```python 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 ) ``` 특징 추출 1024개의 채널 값을 그대로 사용할 수도 있지만, 이는 정보량이 너무 많고 계산 비용이 높아집니다. 따라서, 주요 동위원소의 에너지 값을 추출해 23개의 특징으로 줄였습니다. ```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 isotopes_save(filename: str): with open(filename, "w") as f_out: json.dump(isotopes, f_out) def get_features(spectrum: Spectrum, isotopes: List) -> np.array: energies = [energy for _, energy in isotopes] data = [spectrum.counts[spectrum.energy_to_channel(energy)] for energy in energies] return np.array(data) ``` 모델 훈련 모든 데이터를 하나의 데이터셋으로 결합하고, XGBoost 모델을 사용하여 훈련합니다. GridSearchCV를 통해 최적의 모델 매개변수를 찾습니다. ```python import glob from sklearn.preprocessing import LabelEncoder from xgboost import XGBClassifier from sklearn.model_selection import GridSearchCV all_files = [ ("아메리슘", glob.glob("../data/train/americium.json")), ("라듐", glob.glob("../data/train/radium.json")), ("토륨", glob.glob("../data/train/thorium.json")), ("우라늄 유리", glob.glob("../data/train/uraniumGlass.json")), ("우라늄 유약", glob.glob("../data/train/uraniumGlaze.json")), ("우라니나이트", glob.glob("../data/train/uraninite.json")), ("배경", 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: logging.debug(f"Processing {filename}...") sp = normalize(load_spectrum_json(filename)) for _ in range(augmentation): sp_out = add_noise(sp) x.append(get_features(sp_out, isotopes)) y.append(name) return np.array(x), np.array(y) X_train, y_train = prepare_data(augmentation=10) le = LabelEncoder() le.fit(y_train) y_train = le.transform(y_train) logging.debug("X_train:", X_train.shape) > X_train: (1900, 23) logging.debug("y_train:", y_train.shape) > y_train: (1900,) 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) logging.debug("best_score:", clf.best_score_) > best_score: 0.99474 logging.debug("best_params:", clf.best_params_) > best_params: {'learning_rate': 1.0, 'max_depth': 1, 'n_estimators': 9} isotopes_save("../models/V1/isotopes.json") bst.save_model("../models/V1/XGBClassifier.json") np.save("../models/V1/LabelEncoder.npy", le.classes_) ``` 테스트 모델의 성능을 확인하기 위해, 모델이 미리 보지 못한 데이터를 사용하여 테스트합니다. Radiacode 안드로이드 앱을 통해 XML 파일을 수집했습니다. ```python import xmltodict def load_spectrum_xml(file_path: str) -> Spectrum: with open(file_path) as f_in: doc = xmltodict.parse(f_in.read()) result = doc["ResultDataFile"]["ResultDataList"]["ResultData"] spectrum = result["EnergySpectrum"] cal = spectrum["EnergyCalibration"]["Coefficients"]["Coefficient"] a0, a1, a2 = float(cal[0]), float(cal[1]), float(cal[2]) duration = int(spectrum["MeasurementTime"]) data = spectrum["Spectrum"]["DataPoint"] return Spectrum( duration=duration, a0=a0, a1=a1, a2=a2, counts=[int(x) for x in data], ) test_data = [ ["../data/test/background1.xml", "../data/test/background2.xml"], ["../data/test/thorium1.xml", "../data/test/thorium2.xml"], ["../data/test/uraniumGlass1.xml", "../data/test/uraniumGlass2.xml"], ... ] for group in test_data: data = [] for filename in group: spectrum = load_spectrum_xml(filename) features = get_features(normalize(spectrum), isotopes) data.append(features) X_test = np.array(data) preds = bst.predict(X_test) preds = le.inverse_transform(preds) logging.debug(preds) ``` 모든 결과가 정확하게 나왔습니다. 실시간 테스트 마지막으로, Radiacode 검출기를 사용하여 모델을 실시간으로 테스트해 보겠습니다. ```python from radiacode import RadiaCode, RealTimeData, Spectrum import time import logging le = LabelEncoder() le.classes_ = np.load("../models/V1/LabelEncoder.npy") isotopes = isotopes_load("../models/V1/isotopes.json") bst = XGBClassifier() bst.load_model("../models/V1/XGBClassifier.json") def read_spectrum(rc: RadiaCode): spectrum: Spectrum = rc.spectrum() logging.debug(f"Spectrum: {spectrum.duration} collection time") result = predict_spectrum(spectrum) logging.debug(f"Predict: {result}") def predict_spectrum(sp: Spectrum) -> str: features = get_features(normalize(sp), isotopes) preds = bst.predict([features]) return le.inverse_transform(preds)[0] def read_cps(rc: RadiaCode): data = rc.data_buf() for record in data: if isinstance(record, RealTimeData): logging.debug(f"CPS: {record.count_rate:.2f}") logging.basicConfig(level=logging.DEBUG, format="[%(asctime)-15s] %(message)s", datefmt="%Y-%m-%d %H:%M:%S") rc = RadiaCode() logging.debug(f"ML model loaded") fw_version = rc.fw_version() logging.debug(f"Device connected:, firmware {fw_version[1]}") rc.spectrum_reset() while True: for _ in range(12): read_cps(rc) time.sleep(5.0) read_spectrum(rc) ``` 실제로 방사능 물체를 Radiacode 검출기 근처에 놓고 코드를 실행하면, 모델의 예측이 정확하게 이루어집니다. 결론 이 글에서는 방사능 원소를 예측하는 머신 러닝 모델을 만드는 과정을 설명했습니다. 여러 가지 방사능 샘플을 활용해 모델을 테스트했습니다. 향후 개선 방향으로는 더 많은 데이터 샘플과 동위원소, 더 많은 특징 추가, 다른 모델 유형 테스트 등을 들 수 있습니다. Radiacode 검출기를 사용해 다양한 실험을 할 수 있으며, 이 장비는 저렴한 가격으로 구매 가능합니다. 실제 하드웨어가 없는 독자들은 Kaggle에서 제공되는 데이터를 활용할 수 있습니다. 전체 소스 코드는 Patreon 페이지에서 제공되며, 이 지원은 미래의 테스트를 위한 장비 구매에 도움이 됩니다. LinkedIn을 통해 연결을 맺을 수도 있습니다. 감사합니다.