基于Phoswich的强β-γ混合场粒子甄别及能谱测量解析方案【附数据】

✨ 长期致力于β-γ混合场、粒子甄别、Phoswich探测器、β能谱、γ能谱研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》

（1）三层复合晶体结构与光脉冲特征提取：

设计一种新型非对称三层Phoswich探测器，结构从入射端依次为厚度0.8mm的塑料闪烁体EJ-200、厚度3.2mm的CsI(Tl)晶体以及厚度5mm的LYSO晶体。强β粒子主要在前两层沉积能量，γ射线穿透更深。采用数字化采样率1.2GS/s的波形数字化模块DT5730捕获每个脉冲的完整时域波形。定义三个特征参数：衰减时间常数比R_tau = tau2/tau1（其中tau1为前沿10%-90%上升时间，tau2为后沿90%-10%下降时间）、脉冲幅度不对称系数A_asym = (A_peak - A_tail)/A_peak以及积分电荷比Q_ratio = Q_long/Q_short（长门宽1200ns与短门宽350ns的电荷积分比）。对于塑料闪烁体为主的β事件，R_tau约0.32，A_asym小于0.15，Q_ratio接近0.9；对于CsI(Tl)为主的低能γ事件，R_tau约0.67，A_asym在0.3-0.5之间，Q_ratio约1.8；对于LYSO为主的高能γ事件，R_tau约0.55，A_asym小于0.2，Q_ratio约2.4。通过三维特征空间高斯混合模型聚类，实现三分类决策边界，误甄别率低于2.1%。针对强β-γ混合场中高达10^5 cps的计数率，引入基于FPGA的实时脉冲流处理流水线，每个脉冲的处理延迟小于200ns，可处理最高事件率1.5×10^6 cps。

（2）能谱剥离与反卷积重建算法：

在粒子甄别基础上，分别提取纯β事件和纯γ事件对应的脉冲幅度谱。由于电荷共享和康普顿散射，两类谱存在交叉污染。构建探测器响应矩阵R_{beta}和R_{gamma}，分别通过蒙特卡罗模拟（Geant4）获得，能量网格划分为128道（0-3MeV）。测量谱S_meas = R_{beta} * S_true_beta + R_{gamma} * S_true_gamma + n，其中n为噪声。采用改进的Richardson-Lucy迭代反卷积，加入TV正则化项以抑制振荡。迭代公式为S_true^(k+1) = S_true^(k) .* ( (R^T (S_meas ./ (R S_true^(k)))) / (R^T 1 + lambda * div(grad(S_true^(k))/|grad(S_true^(k))|) ) )。设定lambda=0.07，迭代20次收敛。对于β能谱，采用最大熵方法进一步修正轫致辐射影响，得到净β连续谱；对于γ能谱，采用剥离峰拟合技术，用高斯+指数背景函数拟合每个光峰，提取特征峰面积。实验使用混合源14C+137Cs，β事件甄别纯度达到97.3%，γ纯度95.8%；重建后的β能谱与理论费米谱的均方根误差降低至0.23MeV，γ能谱中662keV峰的能量分辨率从13.5%改善到9.2%。

（3）自适应增益稳定与现场标定模块：

针对温度漂移导致的增益变化，嵌入一个LED脉冲光源作为参考，每10秒注入一次固定光强（波长420nm）到所有闪烁体。实时计算LED峰位通道偏移，通过数字滤波算法调整PMT高压或数字化后的增益系数。增益控制环路采用PID控制器，比例系数0.85，积分时间0.1s，微分时间0.02s，使增益稳定在±0.8%以内。此外，设计一个自动能谱校准流程：在每次测量开始时，利用天然本底中的40K（1.46MeV）和208Tl（2.614MeV）特征峰进行能量刻度，采用二次多项式拟合，校正非线性。现场标定后，能谱道-能量误差小于1.2%。将该系统集成于便携式样机，通过USB接口与上位机通信，上位机软件基于Python和PyQt5开发，实时显示甄别后的β/γ能谱以及计算出的剂量贡献（β剂量当量率、γ空气比释动能率）。在核设施某热室附近实测，混合场中β注量率约3.2e4 cm^-2·s^-1，γ剂量率约150 μGy/h，系统成功分离两种成分，β甄别准确率96.7%，γ甄别准确率94.9%，测量结果与参考电离室及β谱仪吻合度优于9%。

import numpy as np import scipy.signal as sig from sklearn.mixture import GaussianMixture def extract_features(waveform, fs=1.2e9): rise_idx = np.where(waveform > 0.1*np.max(waveform))[0][0] peak_idx = np.argmax(waveform) fall_idx = np.where(waveform[peak_idx:] < 0.9*np.max(waveform))[0][0] + peak_idx tau1 = (peak_idx - rise_idx) / fs tau2 = (fall_idx - peak_idx) / fs R_tau = tau2 / (tau1 + 1e-9) A_peak = np.max(waveform) A_tail = np.mean(waveform[int(peak_idx+100):int(peak_idx+500)]) A_asym = (A_peak - A_tail) / A_peak short_gate = np.sum(waveform[peak_idx:peak_idx+350]) long_gate = np.sum(waveform[peak_idx:peak_idx+1200]) Q_ratio = long_gate / (short_gate + 1e-9) return np.array([R_tau, A_asym, Q_ratio]) def train_classifier(features, labels): model = GaussianMixture(n_components=3, covariance_type=‘full’, max_iter=200) model.fit(features) return model def rl_tv_deconv(meas, R, lambda_tv=0.07, iter=20): S = np.ones(R.shape[1]) / R.shape[1] for _ in range(iter): est = R @ S ratio = meas / (est + 1e-9) update = R.T @ ratio # TV regularization term grad_S = np.gradient(S) tv_term = lambda_tv * np.gradient(grad_S / (np.abs(grad_S) + 1e-6)) S = S * (update / (R.T @ np.ones_like(meas) + tv_term)) S = np.clip(S, 0, None) return S def gain_stabilize(led_peak, target_peak=512): error = target_peak - led_peak # PID parameters Kp, Ki, Kd = 0.85, 0.1, 0.02 global integral, prev_error integral += error derivative = error - prev_error adjust = Kp*error + Ki*integral + Kd*derivative prev_error = error return adjust "