CN103853929B - 一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统及方法 - Google Patents

Abstract

本发明涉及一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析过程及方法，解析过程包括仪器谱探测、建立探测器几何模型、模拟探测器响应函数、响应函数特征参数提取、蒙卡响应矩阵生成、反演解析,依据仪器谱形成的物理过程，建立探测器几何模型，运用蒙卡方法模拟NaI(Tl)闪烁探测器对γ光子的响应函数，确定响应函数的特征参数，并通过插值算法在放射源与γ谱之间构建蒙卡响应矩阵，结合Gold或Boosted-Gold算法，实现在该响应矩阵下反演解析其它被测样品γ仪器谱。应用本发明解析方法省去了谱平滑、寻峰、重峰分解等复杂处理过程，解析结果是待测谱线在该响应矩阵下接近于理论物理谱线的解，该方法对谱线解析的能力提高了。

Description

一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统及方法

技术领域

本发明涉及γ仪器谱解析技术领域，具体涉及一种基于蒙卡(MonteCarlo)响应矩阵的低分辨率γ能谱反演解析系统及方法。

背景技术

在国际原子能机构(IAEA)应用需求的推动下，放射性能谱测量仪已成为核设施等各种核辐射场所、核电站、防化部队、核应急与核查、反核恐怖必不可少的装备。针对机场、港口、海关、车站等大型公共场所，用于环境辐射检测、反核恐怖安检、辐射源清理及天然放射性核素信息探测等核技术应用场合。因此，研究γ能谱解析方法对满足谱仪市场需求、发展核探测与能谱测量技术具有广阔的应用情景和社会现实意义。

由于采用γ能谱仪获得的γ能谱分布与入射到γ探测器之前的γ射线原始谱分布是不同的。通常，把γ射线入射到探测器之前的原始能谱称为γ射线谱，把由γ能谱仪测得的γ能谱称为γ仪器谱，而γ能谱解析正是针对γ仪器谱进行解析的。对于γ射线的能谱测量而言，γ仪器谱是复杂的γ谱。这是因为：首先，由于被测对象本身是多种放射性核素的混合样品，样品放出的γ射线谱是复杂的；其次，γ能谱测量系统受能量分辨本领的限制，尤其是受γ射线探测器的本征能量分辨本领的限制；再次，γ能谱测量系统的环境物体对γ射线的散射本底。相对于半导体探测器而言，NaI(Tl)闪烁探测器具有探测效率高、价格低廉等优势被广泛应用。但由于NaI(Tl)闪烁探测器的能量分辨率有限，使得能量相近的仪器谱峰相互重叠，导致寻峰困难；并且γ光子在NaI(Tl)晶体中产生康谱顿散射使能谱叠加了大量的低能成分，增加了低能区的γ射线总量，造成低能区的谱峰边界模糊，特别在高本底环境下，核素识别率较低甚至错判；进而在对核素种类较多、能谱较复杂样品进行解析时，相应的谱处理算法复杂度也显着增加。对于复杂γ仪器谱的解析，其传统关键技术主要包括谱平滑、寻峰、峰边界确定、本底扣除、重峰分解、净峰面积求取及其活度计算等一系列正演解析过程与方法，由于正演解析过程及方法繁琐，并未考虑探测器与谱仪是否匹配、能谱中多特征参数的提取及放射源与γ能谱之间响应矩阵的构建等自上而下的关联问题，使其放射性核素定性定量分析结果与真实值存在较大差异。

发明内容

本发明的目的是提供了一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统及方法，主要针对NaI(Tl)闪烁探测器能量分辨率低，正演解析过程及方法繁琐，并未考虑探测器与谱仪是否匹配、能谱中多特征参数的提取及放射源与γ能谱之间响应矩阵的构建等自上而下的关联问题，进而实现准确的放射性核素定性定量分析的目标。

为了达到上述目标，本发明提供了基于蒙卡(MonteCarlo)响应矩阵的低分辨率γ能谱反演解析系统，依据核素衰变释放γ光子形成仪器谱的物理过程，建立探测器的几何模型，运用蒙特卡洛方法模拟NaI(Tl)闪烁探测器对γ光子的响应函数，确定探测器响应函数中的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，通过综合仿真和插值算法在放射源与γ能谱之间构建一个蒙卡响应矩阵，结合Gold算法与改进的BoostedGold算法，实现在该响应矩阵下反演解析其它被测样品的γ仪器谱；其特征是：γ仪器谱的解析过程包括仪器谱探测模块、探测器几何模型模块、模拟探测器响应函数模块、响应函数特征参数提取模块、蒙卡响应矩阵生成模块和反演解析模块。

所述仪器谱探测模块，用于将输入的待测样品经NaI(Tl)γ探测谱仪探测后转换为γ仪器谱数据，确定探测器的几何参数，并将所获得的γ仪器谱数据及探测器几何参数分为两路输出，一路输出连接探测器几何模型模块，一路输出连接反演解析模块。

所述探测器几何模型模块，用于根据仪器谱探测模块中的样品谱中的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，并用另一输入标准源的探测器响应函数中全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数予以比对和校正，进而确定探测器的几何模型，使构建的探测器几何模型与被测样品谱数据的测量条件和参数一致。同时，将几何模型参数输出连接至模拟探测器响应函数模块。

所述模拟探测器响应函数模块，用探测器响应函数来反映出全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数和能量的依赖关系，而解析结果的准确性依赖于向解析算法提供的响应矩阵。根据构建的探测器几何模型，采用蒙特卡洛模拟软件模拟探测器响应函数，输出连接至蒙卡响应矩阵生成模块。

所述响应函数特征参数提取模块，用于根据标准源的γ谱来提取该探测器响应函数的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，提取响应函数的特征参数后分两路输出：一路输出连接到探测器几何模型模块，另一路输出连接到蒙卡响应矩阵生成模块。

所述蒙卡响应矩阵生成模块，用于在放射源与γ仪器谱之间构建一个蒙卡响应矩阵，根据模拟探测器响应函数模块输入的能谱数据，通过内插值算法求出其它的响应函数，同时根据响应函数特征参数提取模块输入的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，生成蒙卡响应矩阵，输出连接到反演解析模块。

所述反演解析模块，用于根据蒙卡响应矩阵生成模块输入确定的响应矩阵下，通过反卷积方法解析仪器谱探测模块输入的被测样品γ仪器谱。输出为放射性核素定性定量分析的结果，其解析的结果为测量能谱在蒙卡响应矩阵下对应的能量点(即某个峰位的能量)或接近于理论上的物理谱的解。

本发明所述一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析方法，其特征步骤如下：

步骤1，采用低分辨率NaI(Tl)闪烁探测γ能谱仪，在测量条件和参数(包括探测器晶体直径Φ、探测器的分辨率η、探-源距S₀、源-样距S₁、铅室厚度H等)确定的情况下通过探测待测混合样品，得到待测混合样品的γ仪器谱

步骤2，采用NaI(Tl)闪烁探测γ能谱仪，在与待测样品相同的测量条件下，测量标准源¹³⁷Cs和⁶⁰Co的γ能谱数据，并得到标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)；

步骤3，从步骤2预先建立的标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)与步骤1被测样品谱的特征参数(E₀₁,E₀₂,...E_0m)中选择既定的映射关系比对、校正，进而根据步骤1中实际测量条件和参数构建蒙卡模拟NaI(Tl)闪烁探测器的几何模型M；

步骤4，将步骤3得到的探测器几何模型M，采用MCNP蒙卡模拟软件，通过跟踪单个或多个光子的运动轨迹与反应过程，来实现对响应函数的模拟，得到光子在晶体探测器中的响应函数h(x)；

步骤5，将步骤4所得到模拟光子在晶体探测器中的单个响应函数h(x)，在两条模拟的响应函数之间通过插值算法求出另一响应函数，同时，根据步骤2得到的标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)修正，在50keV～3000keV能量区域内产生所有的响应函数(h₀,h₁,...h_m)，进而得到蒙卡响应矩阵R；

步骤6，根据步骤5获得的蒙卡响应矩阵R并借助解病态逆矩阵的Gold算法来反演解析步骤1测得的混合样品γ仪器谱

本发明的技术方案是构建在放射源与γ谱之间蒙卡响应矩阵，通过该响应矩阵来反演解析低分率NaI(TI)闪烁探测器测量的其它被测样品γ仪器谱。蒙卡响应矩阵解析方法省去了谱平滑、谱寻峰、本底扣除、重峰分解等一系列复杂的传统处理过程，解析得的结果是测量谱在该响应矩阵下对应的能量点或接近于理论上的物理谱的解，且该方法对能谱解析的能力提高了，为谱仪选取合适的探测器及研究能谱数据解析算法提供了理论依据。

附图说明

图1为本发明实施例的系统结构框图。

图2为本发明实施例的方法流程图。

具体实施方式

以下结合附图说明对本发明的实施例作进一步详细描述，但本实施例并不用于限制本发明，凡是采用本发明的相似结构及其相似变化，均应列入本发明的保护范围。

参见图1，本发明实施例提供的基于蒙卡(MonteCarlo)响应矩阵的低分辨率γ能谱反演解析系统，待测样品的γ仪器谱的解析过程包括仪器谱探测模块1、探测器几何模型模块2、模拟探测器响应函数模块3、响应函数特征参数提取模块4、蒙卡响应矩阵生成模块5、反演解析模块6。具体实施时可采用软件固化技术实现各个模块。

所述仪器谱探测模块1，用于将输入的待测样品经NaI(Tl)γ探测谱仪探测后转换为γ仪器谱数据，确定探测器的几何参数，并将所获得的γ仪器谱数据及探测器几何参数分为两路输出，一路输出连接探测器几何模型模块2，一路输出连接反演解析模块6；

所述探测器几何模型模块2，用于根据仪器谱探测模块1中的样品谱中的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，并用另一输入标准源的探测器响应函数中全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数予以比对和校正，进而确定探测器的几何模型，使构建的探测器几何模型与被测样品谱数据的测量条件和参数一致。同时，将几何模型参数输出连接至模拟探测器响应函数模块3。具体实施时，根据建立的几何模型对NaI(Tl)闪烁探测器点源¹³⁷Cs进行模拟，在蒙卡模拟软件MCNP中，不直接对几何体进行描述，而是通过对围成该几何体的面进行描述来时实现几何体或几何体区域的定义。例如：MCNP输入文件的结构中的几何区域定义，曲面的几何位置，粒子在区域中的传输，包括源的位置、能量、能谱、方向等等，由于γ射线的吸收与散射，标准源谱与物理谱存在很大的差异，如果探测器响应函数不给予修正，那么从γ谱分析中获得的也是不准确的物理数据。因此，根据样品谱的测量条件和标准源的测量参数进行比对校正，进而得到了蒙卡模拟NaI(Tl)闪烁探测器几何模型。为了和实验测量条件尽量接近，模拟计算模型中还包含闪烁晶体、反射层，包裹层以及光导(SiO₂)。

所述模拟探测器响应函数模块3，用探测器响应函数来反映出全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数和能量的依赖关系，而解析结果的准确性依赖于向解析算法提供的响应矩阵。根据构建的探测器几何模型，采用蒙特卡洛模拟软件模拟探测器响应函数，输出连接至蒙卡响应矩阵生成模块5。具体实施时，根据γ能谱对应的能量沉积是跟踪大量γ源光子从产生到进入晶体探测器再到光子被吸收或逸出探测器的整个光子-电子输运过程。能量为E的光子从源发出，入射到探测器上，并在探测器内发生各种作用沉积能量。例如：从源S处发出能量为E₀的光子，入射到NaI(Tl)探测器上并在探测器晶体内发生各种相互作用，能量沉积，激发荧光，从而使探测器内沉积的能量激发出荧光引起幅度响应脉冲。因此，响应函数是描述入射光子能量与其引起的脉冲幅度之间关系的随机关系函数，进而反映出全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数和能量的依赖关系。

所述响应函数特征参数提取模块4，用于根据标准源的γ谱来提取该探测器响应函数的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，提取响应函数的特征参数后分两路输出：一路输出连接到探测器几何模型模块2，另一路输出连接到蒙卡响应矩阵生成模块5；实施例的计算过程为，在用蒙卡MCNP软件模拟时，沉积能量的计算是把每个源粒子事件在探测器中发生相互作用的能量沉积累积以后得到。考虑到谱仪系统分辨率，沉积能量的高斯展宽，其谱仪半宽度FWHM和射线能量E_d间的关系可以写成如下形式：

$FWHM\left(E_d\right)=a+b\sqrt{E_d+{\mathrm{cE}}_d^2}---\left(1\right)$

式(1)中a，b，c是分辨刻度系数，取a＝0.001，b＝0.05086，c＝0.030486。记录能量。

E′＝E_d+σ·x(2)

式(2)中，x由标准正态分布抽样得到。

在计算谱仪的探测效率时，模拟时需要记录全谱的总计数n：包括探测器中的γ光子总数N及被探测器记录的光子总数，利用全能峰下的计数n_p与全谱的总计数n之比可计算出峰总比R[E]。

所述蒙卡响应矩阵生成模块5，用于在放射源与γ仪器谱之间构建一个蒙卡响应矩阵，根据模拟探测器响应函数模块3输入的谱数据，通过内插值算法求出其它的响应函数，同时根据响应函数特征参数提取模块4输入的全能峰，康普顿边缘，康普顿平台，反散射峰，逃逸峰等特征参数，生成蒙卡响应矩阵，输出连接到反演解析模块6；具体实施时，是基于反卷积迭代的解析方法，谱解析是对线性系统方程求解：

y＝R·x(3)

式(3)中，y表示被测量的能谱，R表示响应矩阵，x表示分解后生成的能谱。假设整个能谱段的光电峰形状都是相同的。通过对全能峰响应函数中每个相邻列的转置构建响应矩阵R。然后由反卷积算法可获得γ谱的反卷结果。

所述反演解析模块6，用于根据蒙卡响应矩阵生成模块5输入确定的响应矩阵下，通过反卷积方法解析仪器谱探测模块1输入的被测样品γ仪器谱。输出为放射性核素定性定量分析的结果，其解析的结果为测量谱在蒙卡响应矩阵下对应的能量点(即某个峰位的能量)或接近于理论上的物理谱的解。具体实施时，通过在放射源与γ谱之间构建蒙卡响应矩阵R去反演解析被测的样品谱，但由于直接解蒙卡响应矩阵方程反卷积结果受干扰或噪声影响很大，不可能通过解病态矩阵方程得到真值，因此，通过采用Gold和Boosted-Gold非线性迭代逼近稳定点的方法，并把收敛到的稳定值当成矩阵的解。

参见图2，本发明实施例提供的基于蒙卡(MonteCarlo)响应矩阵的低分辨率γ能谱反演解析方法，可采用计算机软件技术手段自动进行流程，具体包括以下步骤：

步骤1，采用低分辨率NaI(Tl)闪烁探测γ能谱仪，在测量条件和参数(包括探测器晶体直径Φ、探测器的分辨率η、探-源距S₀、源-样距S₁、铅室厚度H等)确定的情况下通过探测待测混合样品，得到待测混合样品的γ仪器谱实施例中，将待测样品铀、钍、钾分别300g放入圆柱型样品盒放置探测器前端，测量置于铅室中，铅室厚度为50mm，采用的探测器为NaI(Tl)闪烁探测器，直径Φ75mm×75mm，自研制的1024道DSP数字多道谱仪，探测器的能量分辨率7.9％，积分非线性<0.1％，探测器与被测样品距离为50mm，能量范围50～3000keV，计数率>20k，自带稳谱功能，放射性测量采用相对测量方法，得到了1024道的被测样品γ仪器谱

步骤2，采用NaI(Tl)闪烁探测1024道γ能谱仪，在与待测样品相同的测量条件下，测量标准源¹³⁷Cs和⁶⁰Co的γ能谱数据，并得到标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)。实施例中，根据核素放射出γ光子形成仪器谱的物理过程，采用用两种标准源的γ谱为¹³⁷Cs-661keV，⁶⁰Co-1173keV，1332keV，通过对比仿真的响应函数和测量标准源的响应函数，来确定NaI(Tl)闪烁探测器仿真最佳的几何角度，并获得¹³⁷Cs-661keV标准源γ能谱响应函数的特征参数，例如:E₁₁-全能峰(661keV)，E₁₂-康普顿边缘(480keV)，E₁₃-康普顿平台(200keV～480keV)，E₁₄-反散射峰(185keV)，E₁₅-逃逸峰(511keV和1022keV)等特征参数。对于能量线性关系的检验，是用已知能量的一组标准源，测量相应全能峰处的脉冲幅度，建立γ射线能量与对应峰位的关系曲线，即刻度曲线，是一条不通过原点的直线，其公式：

E(X_p)＝GX_p+E₀(4)

式(4)中，X_p为全能峰峰位，E₀为直线截距，G为增益(单位脉冲幅度对应的能量)。

步骤3，从步骤2预先建立的标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)与步骤1被测样品谱的特征参数(E₀₁,E₀₂,...E_0m)中选择既定的映射关系比对，进而根据步骤1中实际测量条件和参数构建蒙卡模拟NaI(Tl)闪烁探测器的几何模型M。实施例中，为了和实验测量条件尽量接近，探测器的结构由衬底面直径75mm，高75mm的NaI(Tl)闪烁圆柱体，在闪烁圆柱体内包一层厚3mm的铝壳，后面放有底直径为75mm，高为20mm的耦合光电倍增管等部分组成。全部放置于一个内半径90mm，外半径140mm内、外高分别为300mm和360mm的铅屏蔽层中，模拟计算模型包含闪烁晶体、反射层，包裹层以及光导(SiO₂)。晶体尺寸为Φ75mm×75mm，晶体前部A_l壳厚3mm，外侧面A_l壳厚2mm，MgO₂反射层厚0.5mm，后部为厚2mm的SiO₂光学玻璃。

步骤4，将步骤3得到的探测器几何模型M，采用MCNP蒙卡模拟软件，通过跟踪单个或多个光子的运动轨迹与反应过程，来实现对响应函数的模拟，得到光子在晶体探测器中的响应函数h(x)。响应函数是描述入射光子能量与其引起的脉冲幅度之间关系的随机关系函数，仿真模拟的主要目标是产生对应几何模型尽可能接近实际探测器的响应函数，实施例中，根据建立的几何模型对NaI(Tl)闪烁探测器点源¹³⁷Cs进行模拟。在MCNP中，不直接对几何体进行描述，而是通过对围成该几何体的面进行描述来时实现几何体或几何体区域的定义。

NaI(Tl)探测器点源¹³⁷Cs模拟

将直接输出由蒙特卡洛模拟后获得的γ光子在晶体探测器中沉积的能量，经过统计得到每种能量与其对应的光子数，并对每个点按高斯分布展宽后即可获得¹³⁷Cs的γ能谱响应函数h(x)。

步骤5，将步骤4所得到模拟光子在晶体探测器中的单个响应函数h(x)，在两条模拟的响应函数之间通过插值算法求出另一响应函数，同时，根据步骤2得到的标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)修正，在50keV～3000keV能量区域内产生所有的响应函数(h₀,h₁,...h_m)，进而得到蒙卡响应矩阵R。实施例中，通过蒙卡模拟软件MCNP，利用NaI(Tl)闪烁探测器几何模型，对单能¹³⁷Cs在50keV-3000keV能量范围内的入射γ射线进行模拟，由于响应函数是可以通过几种单能γ谱获得的，因此必须通过对几种单能γ谱的响应函数进行仿真，才能得到多能段谱，但响应函数的仿真很耗时，因此在实施例中，采用每隔100keV获得一个响应函数，两个模拟的响应函数间通过内插值算法求出另一个响应函数。响应函数的插值算法描述如下：

首先，假设在不同能量段有两个响应函数分别是h₁(E₁₁,e),h₂(E₁₂,e)，将假定插值段分为n部分。然后，对于每个插值响应函数h_i(E_1i,e)的点可分类计算得到。如下：

$e_i\left(k\right)=\frac{E_{1i}-E_{12}}{E_{11}-E_{12}}\&lsqb;e_1\left(k\right)-e_2\left(k\right)\&rsqb;+e_2\left(k\right),k\&Element;<1,n>---\left(5\right)$

对于每部分k∈{1,n}和每个e_i∈{e_i(k-1),e_i(k)}，找出匹配的一对点(e₁(j_k),e₂(j_k))：

$e_1\left(j_k\right)=\frac{e_i\left(j_k\right)-e_i(k-1)}{e_i\left(k\right)-e_i(k-1)}\&lsqb;e_1\left(k\right)-e_1(k-1)\&rsqb;+e_1(k-1)---\left(6\right)$

$e_2\left(j_k\right)=\frac{e_i\left(j_k\right)-e_i(k-1)}{e_i\left(k\right)-e_i(k-1)}\&lsqb;e_2\left(k\right)-e_2(k-1)\&rsqb;+e_2(k-1)---\left(7\right)$

从响应函数h₁,h₂中，读出其值为：

y_i(j_k)＝h₁(E₁₁,e₁(j_k)),y₂(j_k)＝R₂(E₁₂,e₂(j_k))(8)

然后，响应矩阵h_i中对应能量e_i的值由插值给出如下：

$y_i\left(j_k\right)=\frac{E_{1i}-E_{12}}{E_{11}-E_{12}}\&lsqb;y_1\left(j_k\right)-y_2\left(j_k\right)\&rsqb;+y_2\left(j_k\right)---\left(9\right)$

用上述算法在给定的区域内可生产所有的响应函数(h₀,h₁,...h_m)。通过填充所有的间隔区域，按列分布可以获得整个响应矩阵R。即：

在探测器中，由于γ射线的吸收与散射，标准源谱与物理谱存在很大的差异，因此，如果探测器响应函数不给予修正，那么从γ谱分析中获得的也是不准确的物理数据。实施例中，按步骤2特征参数定义次序将响应函数分成了4部分：(0,E₁₄)，(E₁₄,E₁₂)，(E₁₂,E₁₁-10keV)和(E₁₁-10keV,E₁₁+10keV)来修正响应函数，其中E₁₄＝E₁₁/(1+2E₁₁/m_ec²)，E₁₂＝E₁₁-E₁₄，E₁₁分别是反散射峰，康普顿边缘和全能峰的能量。当能量低于E＝m_ec²/2时，康普顿边缘和反散射峰重叠，区域选择变成为(0,E₁₄)，(E₁₄,E₁₂)，(E₁₂,E₁₁-10keV)和(E₁₁-10keV,E₁₁+10keV)，此处E₁₂＜E₁₁。在200～300keV范围内混叠了康普顿边缘和反散射峰，为了达到更好的分辨率和插值方法的准确性，通过仿真模拟获得响应函数的间隔缩小到50keV。对于能量超过1.02MeV电子对效应凸显，结果导致在γ能谱上出现逃逸峰。因此，当能量超过1.02MeV时进行响应函数的内插值时必须要特别注意。逃逸峰被单独分开进行插值，最终为每个插值响应函数中将逃逸峰加到产生的准确位置(E-511keV和E-1022keV)上。

步骤6，根据步骤5获得的蒙卡响应矩阵R并借助解病态逆矩阵的Gold算法来反演解析步骤1测得的混合样品γ仪器谱连续线性系统中卷积公式常表示为：y(t)＝h(t)*x(t)，其中x(t)为输入函数，h(t)为系统函数，y(t)为输出函数，三个量中知道任意两个就能求出第三个。实施例中，首先，将γ仪器谱的成谱可表示为：

$\widetilde{y_i\left(n\right)}=h\left(n\right)*x\left(n\right)---\left(11\right)$

谱的解析过程变成了由所测得的仪器谱和预知的γ谱探测器响应h(n)求解输入x(n)的过程。将式(11)用矩阵方程表示：

式(12)中向量为n道的仪器谱，向量为m道，与仪器谱不同，这里每道对应物理谱的一个能量点，即入射γ射线的能量；R为响应矩阵，列向量R_～,j表示第i能量点对应的响应，同样是一条n道的能谱，i＝1…m。

其次，在实施例中反算法的实现。由式可知但不论是方阵求逆的苛刻条件、矩阵的病态性还是线性方程组的复杂性，都决定了直接求取x是不可行性的。根据仪器谱非负定的物理特点，采用迭代逼近稳定点的方法，将收敛到的一个精确稳定点作为解。选取Gold迭代反卷积算法和改进型的Boosted-Gold算法作为反卷积谱解析算法。反卷积求解迭代的具体过程如下：

首先，在式(12)左右两边同时乘以R^T得到系数为常对角矩阵的方程：

z＝Ax(13)

式(13)中，A＝R^TR，z＝R^Ty。则第k+1次迭代的结果可以表示为：

x^(k+1)＝x^(k)+μ(z-Ax)(14)

μ为迭代过程中引入的局部变量的松弛因子：

${\mathrm{\&mu;}}_i=\frac{x_i^{\left(k\right)}}{{\mathrm{\&Sigma;}}_{m=1}^N{A}_{im}{x}_m^{\left(k\right)}}---\left(15\right)$

将局部变量松弛因子μ_i代入式(15)得：

$x_i^{(k+1)}=x_i^{\left(k\right)}+\frac{x_i^{\left(k\right)}}{{\mathrm{\&Sigma;}}_{m=1}^N{A}_{im}{x}_m^{\left(k\right)}}\&lsqb;z_i-\underset{m=1}{\overset{N}{\&Sigma;}}{A}_{im}{x}_m^{\left(k\right)}\&rsqb;---\left(16\right)$

将式(16)化简得：

$x_i^{(k+1)}=\frac{z_i}{{\mathrm{\&Sigma;}}_{m=1}^N{A}_{im}{x}_m^{\left(k\right)}}{x}_i^{\left(k\right)},i=1,2,...,N,k=0,1,...,L---\left(17\right)$

式(17)即为Gold算法的最终迭代式，将迭代初值设为:

x⁽⁰⁾＝[1,1,...,1]^T(18)

将式(18)作为初值开始迭代，最终将收敛于稳定点集，即x^(k+1)＝＝x^(k)，所得的稳定点集就是反演结果，即为所测得的样品谱解析的结果。

利用Gold解析收敛到稳定值以后，无论如何增加迭代次数，其解析结果也不会再改变了。继续缩小峰的宽度发现，当解到达稳定状态时需要停止迭代，然后通过某种方式去改变向量x^(L)，以此为新的初值，再重复Gold解谱迭代式(17)。并需要应用非线性增强(加速)函数去改变特解。幂函数被证实可以给出较好的结果。其增强改进的迭代Boosted-Gold解析算法如下所述：

①根据式(18)设置初始值为x⁽⁰⁾；②根据需要设置重复次数R和迭代次数L；③使重复次数值r＝1；④根据式(17)，针对k＝0,1,...,L-1寻找解x^(L)；⑤如果r＝R，终止计算。否则继续下列运算：a)进行加速增强操作，即其中p是加速指数；b)使r＝r+1；c)继续重复第⑤步。

上述的实施例中，运用该反演解析方法对²³⁸U系、²³²Th系、⁴⁰K、²³⁸U-²³²Th混合源等样品谱的解析，解析得的结果是测量谱在该响应矩阵下对应的能量点或接近于理论上的物理谱的解，且该方法对能谱解析的能力提高了。结果表明蒙卡响应矩阵解析方法减小了统计涨落的影响，且可省去谱平滑、谱寻峰、本底扣除、重峰分解等一系列复杂的处理过程，其成果为谱仪选取合适的探测器及研究能谱数据解析算法提供了理论依据。

Claims (8)

1.一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，依据核素衰变释放γ光子形成仪器谱的物理过程，建立探测器的几何模型，运用蒙特卡洛方法模拟NaI(Tl)闪烁探测器对γ光子的响应函数，确定探测器响应函数中的特征参数，通过综合仿真和插值算法在放射源与γ能谱之间构建一个蒙卡响应矩阵，结合Gold算法与改进的BoostedGold算法，实现在该响应矩阵下反演解析其它被测样品的γ仪器谱；其特征是：γ仪器谱的解析过程包括仪器谱探测模块、探测器几何模型模块、模拟探测器响应函数模块、响应函数特征参数提取模块、蒙卡响应矩阵生成模块和反演解析模块。

2.根据权利要求1所述的一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，其特征在于：所述仪器谱探测模块，用于将输入的待测样品经NaI(Tl)γ探测谱仪探测后转换为γ仪器谱数据，确定探测器的几何参数，并将所获得的γ仪器谱数据及探测器几何参数分为两路输出，一路输出连接探测器几何模型模块，一路输出连接反演解析模块。

3.根据权利要求1所述的一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，其特征在于：所述探测器几何模型模块，用于根据仪器谱探测模块中的样品谱中的特征参数，并用另一输入标准源的探测器响应函数中特征参数予以比对和校正，进而确定探测器的几何模型，使构建的探测器几何模型与被测样品谱数据的测量条件和参数一致；同时，将几何模型参数输出连接至模拟探测器响应函数模块。

4.根据权利要求1所述的一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，其特征在于：所述模拟探测器响应函数模块，用探测器响应函数来反映出特征参数和能量的依赖关系，而解析结果的准确性依赖于向解析算法提供的响应矩阵；根据构建的探测器几何模型，采用蒙特卡洛模拟软件模拟探测器响应函数，输出连接至蒙卡响应矩阵生成模块。

5.根据权利要求1所述的一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，其特征在于：所述响应函数特征参数提取模块，用于根据标准源的γ谱来提取该探测器响应函数的特征参数，提取响应函数的特征参数后分两路输出：一路输出连接到探测器几何模型模块，另一路输出连接到蒙卡响应矩阵生成模块。

6.根据权利要求1所述的一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，其特征在于：所述蒙卡响应矩阵生成模块，用于在放射源与γ仪器谱之间构建一个蒙卡响应矩阵，根据模拟探测器响应函数模块输入的谱数据，通过内插值算法求出其它的响应函数，同时根据响应函数特征参数提取模块输入的特征参数，生成蒙卡响应矩阵，输出连接到反演解析模块。

7.根据权利要求1所述的一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析系统，其特征在于：所述反演解析模块，用于根据蒙卡响应矩阵生成模块输入确定的响应矩阵下，通过反卷积方法解析仪器谱探测模块输入的被测样品γ仪器谱；输出为放射性核素定性定量分析的结果，其解析的结果为测量谱在蒙卡响应矩阵下对应的能量点或接近于理论上的物理谱的解。

8.一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析方法，其特征步骤如下：

步骤1，采用低分辨率NaI(Tl)闪烁探测γ能谱仪，在测量条件和参数确定的情况下通过探测待测混合样品，得到待测混合样品的γ仪器谱

步骤2，采用NaI(Tl)闪烁探测γ能谱仪，在与待测样品相同的测量条件下，测量标准源¹³⁷Cs或⁶⁰Co的γ能谱数据，并得到标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)；

步骤3，从步骤2预先建立的标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)与步骤1被测样品谱的特征参数(E₀₁,E₀₂,...E_0m)中选择既定的映射关系比对、校正，进而根据步骤1中实际测量条件和参数构建蒙卡模拟NaI(Tl)闪烁探测器的几何模型M；

步骤4，将步骤3得到的探测器几何模型M，采用MCNP蒙卡模拟软件，通过跟踪单个或多个光子的运动轨迹与反应过程，来实现对响应函数的模拟，得到光子在晶体探测器中的响应函数h(x)；

步骤5，将步骤4所得到模拟光子在晶体探测器中的单个响应函数h(x)，在两条模拟的响应函数之间通过插值算法求出另一响应函数，同时，根据步骤2得到的标准源γ能谱响应函数的特征参数(E₁₁,E₁₂,...E_1m)修正，在50keV～3000keV能量区域内产生所有的响应函数(h₀,h₁,...h_m)，进而得到蒙卡响应矩阵R；

步骤6，根据步骤5获得的蒙卡响应矩阵R并借助解病态逆矩阵的Gold算法来反演解析步骤1测得的混合样品γ仪器谱