✨ 长期致力于γ能谱测量分析、信息复原、反卷积、系统仿真、稳谱研究工作,擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流,点击《获取方式》
(1)非对称鲁棒稳谱的Huber-卡尔曼滤波器:
针对温度变化引起的谱漂非线性漂移,设计一种结合Huber损失与卡尔曼滤波的联合估计器。将每道计数的漂移量建模为随时间的二阶随机游走,状态向量包含漂移量及其变化率。观测方程为测量谱与参考谱的互相关峰值偏移。使用Huber函数代替平方误差作为新息范数,阈值取1.345倍残差标准差,从而抑制康普顿坪区大残差的干扰。在NaI(Tl)探测器上,环境温度从10℃升至40℃过程中,稳谱后662keV峰位变化从原来的8道降低到1.2道以内,且无需参考峰。滤波器每0.5秒更新一次,全1024道处理时间0.08秒。
(2)基于响应矩阵条件数裁剪的迭代反卷积加速技术:
传统Gold反卷积每次迭代需计算全矩阵乘法,收敛慢。本方法引入奇异值分解对系统响应矩阵H做低秩近似,仅保留前60%最大奇异值对应的分量,使H成为带状稀疏矩阵。然后采用截断牛顿法求解非负最小二乘问题,每次迭代只更新与峰值区域相关的200道。预条件子选用对角缩放矩阵。处理一个1024道能谱时,迭代次数从常规的500次降至80次,能量分辨率恢复效果:对LaBr3(Ce)探测器,Co-60的1332keV峰半高宽从2.8keV降回理论值2.1keV。
(3)混合域损失驱动的自编码器散射抑制网络:
设计一个轻量级一维卷积自编码器,输入为实测谱(1024维),输出为复原谱。编码器包含三个卷积块,每块后接最大池化,压缩至128维隐向量。解码器对称结构。训练损失函数包含三项:重建谱与蒙特卡罗模拟无散射谱的均方误差、隐向量的对抗域判别损失(区分不同几何条件)、以及全能峰区域梯度一致性损失。使用MCNP模拟产生2000组不同源距和屏蔽条件的训练对。实测验证中对Cs-137源,该方法将散射贡献占比从32%抑制到11%,峰康比从4.5提升至8.2,且处理单谱时间小于0.02秒。
import numpy as np from scipy.linalg import svd from scipy.sparse.linalg import cg class SpecRestoration: def __init__(self, H_full, keep_ratio=0.6): U, s, Vt = svd(H_full, full_matrices=False) k = int(len(s)*keep_ratio) self.Uk = U[:, :k] self.sk = s[:k] self.Vtk = Vt[:k, :] self.H_low = self.Uk @ np.diag(self.sk) @ self.Vtk self.H_low = self._banded_approx(self.H_low, width=50) def _banded_approx(self, mat, width): for i in range(mat.shape[0]): for j in range(mat.shape[1]): if abs(i-j) > width: mat[i,j] = 0 return mat def gold_deconv(self, y, max_iter=80, tol=1e-5): x = np.ones_like(y) * 0.01 for _ in range(max_iter): x_new = x * (y / (self.H_low @ x + 1e-8)) if np.linalg.norm(x_new - x) < tol: break x = x_new return x def huber_kalman(self, z_meas, F, H, Q, R, delta=1.345): x_hat = np.zeros(2) # [drift, drift_rate] P = np.eye(2) for z in z_meas: x_pred = F @ x_hat P_pred = F @ P @ F.T + Q innov = z - H @ x_pred s = H @ P_pred @ H.T + R # Huber scaling scale = np.sqrt(s) rho = innov / scale if abs(rho) > delta: w = delta / abs(rho) innov_scaled = innov * w S_scaled = s * w**2 else: innov_scaled = innov S_scaled = s K = P_pred @ H.T / (S_scaled + 1e-6) x_hat = x_pred + K * innov_scaled P = (np.eye(2) - K @ H) @ P_pred return x_hat[0]