✨ 长期致力于核废物桶、层析γ扫描、图像重建、压缩感知研究工作,擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流,点击《获取方式》
(1)压缩感知IART-TVM算法与极坐标网格划分:
针对单HPGe探测器TGS系统,提出压缩感知IART-TVM图像重建算法。将核废物桶划分为极坐标网格(12圈×72块,共864体素),比直角坐标网格减少约40%体素数量。测量采用稀疏采样:水平位置选4个角度(0°, 90°, 180°, 270°),每个位置24个角度步(步长15°),共96组投影数据,扫描速度提高为原来的9倍。建立投影数据与松弛因子之间的二项式关系:lambda_k = 0.2 + 0.8/(1+exp(-k/10))。代数重建算法ART迭代中,引入全变分最小化TVM作为正则项,形成IART-TVM。在仿真非均匀核废物桶(7种密度介质+Cs-137源)中,IART-TVM重建的线衰减系数均方差为0.028,而ART为0.044;信噪比提升64%。核素活度重建值误差3.21%,定位准确(体素偏差不超过1个网格)。
(2)MLEM-TVM优化算法提升初始图像质量:
将极大似然期望最大化MLEM算法用于生成初始图像,因其保真度高。然后结合TVM迭代细化。MLEM迭代20次后得到初始图像,再运行TVM 30次。比较表明,MLEM-TVM的均方差为0.021,比IART-TVM的0.028降低25%;信噪比2.34倍于纯MLEM。活度重建误差3.01%,比IART-TVM的4.35%更低。在投影数据欠采样(仅48组)条件下,MLEM-TVM仍能将信噪比维持在18dB以上,而ART降为11dB。该算法对噪声鲁棒,在计数统计波动±15%时,重建活度变化<5%。
(3)点源空间效率函数TGS效率刻度优化:
利用MCNP建立TGS探测系统模型,计算507个空间点(11层×13×13)对多能量射线的探测效率。提出点源空间效率函数模型:epsilon(r, E) = (a1*exp(-b1*r) + a2*exp(-b2*r)) * (c0 + c1*E + c2*E^2)。对59keV, 81keV, 662keV, 1173keV, 1332keV拟合,相关系数R^2均>0.99。体素效率取体素中心点效率近似,误差<6%。编写Matlab GUI软件TGS-Rec,集成能谱分析(自动寻峰、净峰面积计算)、效率刻度、图像重建功能。处理实际废物桶时,Cs-137总活度测量值3.143e5 Bq,与真值误差2.71%。
import numpy as np from scipy.sparse.linalg import lsqr from skimage.restoration import denoise_tv_chambolle def iart_tvm_reconstruction(projections, angles, n_voxels=864, n_iter=30): # 压缩感知IART-TVM算法 # projections: 测量投影数据 # angles: 对应角度 x = np.zeros(n_voxels) for k in range(n_iter): # 松弛因子二项式关系 lambda_k = 0.2 + 0.8 / (1 + np.exp(-k/10)) # ART迭代 (简化) for i, (p, theta) in enumerate(zip(projections, angles)): # 前向投影 Ai = np.random.randn(n_voxels) # 投影矩阵行 px = Ai @ x if np.linalg.norm(Ai) > 0: x += lambda_k * (p - px) / (np.linalg.norm(Ai)**2) * Ai # 全变分最小化 x_img = x.reshape(24, 36) # 假设12圈x72块变形成2D x_tv = denoise_tv_chambolle(x_img, weight=0.1) x = x_tv.ravel() return x def mlem_tvm(projections, system_matrix, n_iter_mlem=20, n_iter_tvm=30): # MLEM-TVM 优化算法 x = np.ones(system_matrix.shape[1]) / system_matrix.shape[1] for _ in range(n_iter_mlem): # MLEM更新 forward = system_matrix @ x ratio = projections / (forward + 1e-6) x = x * (system_matrix.T @ ratio) / (system_matrix.T @ np.ones_like(projections)) x = np.maximum(x, 1e-6) for _ in range(n_iter_tvm): x_img = x.reshape(24, 36) x_tv = denoise_tv_chambolle(x_img, weight=0.05) x = x_tv.ravel() return x def point_source_efficiency(distance, energy, coeff): # 点源效率函数 # coeff: [a1, b1, a2, b2, c0, c1, c2] a1,b1,a2,b2,c0,c1,c2 = coeff r_term = a1 * np.exp(-b1 * distance) + a2 * np.exp(-b2 * distance) e_term = c0 + c1 * energy + c2 * energy**2 return r_term * e_term def tgs_gui_spectrum_analysis(spectrum, energy_calib): # 能谱分析:寻峰和净面积计算 from scipy.signal import find_peaks peaks, _ = find_peaks(spectrum, height=100, distance=10) net_areas = [] for pk in peaks: left = max(0, pk-5) right = min(len(spectrum), pk+6) background = np.mean(np.concatenate([spectrum[left:pk-2], spectrum[pk+2:right]])) net_area = np.sum(spectrum[pk-2:pk+3]) - background * 5 net_areas.append(net_area) return peaks, net_areas # 示例 if __name__ == '__main__': # 模拟投影 n_proj = 96 n_vox = 864 proj_sim = np.random.randn(n_proj) + 10 rec = iart_tvm_reconstruction(proj_sim, np.linspace(0, 360, n_proj)) print('Reconstruction shape:', rec.shape) # 效率刻度 coeffs = [0.45, 0.08, 0.32, 0.15, 0.02, 0.001, -1e-7] eff = point_source_efficiency(distance=0.5, energy=0.662, coeff=coeffs) print('Efficiency at 50cm, 662keV:', eff)