CN103135125A - 一种过滤x射线参考辐射的能谱分析方法 - Google Patents

Abstract

本发明属于电离辐射计量技术，具体涉及一种过滤X射线参考辐射的能谱分析方法。本发明用HPGe能谱仪测量低空气比释动能率和窄谱两个系列的过滤X射线参考辐射，探索出使用反卷能谱法对测量谱进行了能谱分析的方法，得到了参考辐射场的注量谱和各特征值，并与国标中给出的参考谱和各特征值的推荐值进行了比较，两者符合的较好，验证了已建立的X射线参考辐射的准确性和可靠性。

Description

一种过滤X射线参考辐射的能谱分析方法

技术领域

本发明属于电离辐射计量技术，具体涉及一种过滤X射线参考辐射的能谱分析方法。

背景技术

过滤X射线参考辐射在电离辐射计量中是非常重要的参考辐射。国家标准″GB/T 12162-2000″规定了不同系列的过滤X射线参考辐射的平均能量、能量分辨率、半值层、同质系数等参数，并推荐通过能谱研究证实平均能量和能量分辨率。连续X射线能谱分析一直是核探测领域的难题，早期是使用NaI探测器进行实验测量，用剥谱法等方法进行能谱分析，得到的能谱精度较低；HPGe谱仪发展起来后，国际上开始用HPGe谱仪连续测量X射线能谱，用剥谱法、最小二乘法或反卷能谱法分析测量谱得到的注量谱；而在国内，虽然GB/T 12162-2000中推荐用能谱研究的方法确定过滤X参考辐射的特征参数，但是近几年国内外对相关工作鲜有报道。

发明内容

本发明的目的在于提供一种用HPGe能谱仪测量低空气比释动能率和窄谱两个系列的过滤X射线参考辐射，并对测量谱进行能谱分析的方法，得到参考辐射场的注量谱和各特征值，用以验证已建立的X射线参考辐射的准确性和可靠性。

本发明的技术方案如下：一种过滤X射线参考辐射的能谱分析方法，包括如下步骤：

(1)建立由探测器、数字化谱仪、计算机组成的能谱测量系统，并对能谱测量系统进行能量刻度，测量得到过滤X射线的测量谱；

(2)用蒙特卡罗模拟和数学插值的方法确定能谱测量系统的响应矩阵；

(3)对测量得到的测量谱进行剥谱法和反卷能谱法分析，得到了过滤X射线的注量谱，并与理论计算谱进行比较。

进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，所述的能谱测量系统的探测器为HPGe探测器。

进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，步骤(1)中对能谱测量系统进行能量刻度的方法如下：取一同位素刻度源，接近探测器放置，开始测量其能谱；一段时间后，测量谱中出现明显的峰，将峰中心能量对应道址的能量改为该同位素刻度源的能量；再另取一不同能量的同位素刻度源，用相同的办法进行能量刻度；反复进行几次，挑选的同位素刻度源能量应能够覆盖即将测量连续X射线的能量范围。

更进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，步骤(1)中对能谱测量系统进行能量刻度所选用的同位素刻度源包括：¹³⁷Cs、²⁴¹Am、¹⁵²Eu、¹³³Ba。

进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，步骤(2)中所述的蒙特卡罗模拟方法是首先对HPGe探测器及周围的其他实验装置和环境进行数学建模，然后通过模拟大量由放射源发射并通过HPGe探测器的单个X光子的历史，统计每个X光子在HPGe探测器中沉积的能量，得到服从统计分布的在HPGe探测器中能量沉积的分布。

更进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，在蒙特卡罗模拟计算中使用MCNP5程序提供的F8计数类型卡记录光子在晶体中沉积的能量脉冲的幅度分布，然后选用计数特殊处理卡GEB对脉冲能量幅度分布进行高斯展宽。

进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，步骤(2)中所述的数学插值的方法是采用分段内插的方式，得到能量范围为6keV～305keV的光子在能谱测量系统中的完整的响应矩阵。

进一步，如上所述的过滤X射线参考辐射的能谱分析方法，其中，步骤(3)中从测量谱得到注量谱是通过Gold反卷算法求解线性方程组M＝R×F，

其中，F为注量谱，R为响应矩阵，M是实验测量得到的测量谱。

本发明的有益效果如下：本发明用HPGe能谱仪测量低空气比释动能率和窄谱两个系列的过滤X射线参考辐射，探索出使用反卷能谱法对测量谱进行了能谱分析的方法，得到了参考辐射场的注量谱和各特征值，并与国标中给出的参考谱和各特征值的推荐值进行了比较，两者符合的较好，验证了已建立的X射线参考辐射的准确性和可靠性。

附图说明

图1为本发明的方法流程图；

图2为能谱测量系统的实验装置示意图；

图3为确定探测器晶体尺寸的实验示意图；

图4为¹³³Ba的模拟能谱和实测能谱比较图；

图5为解谱程序流程图；

图6为低空气比释动能率系列解谱得到的注量谱和理论谱的比较示意图；

图7为窄谱系列解谱得到的注量谱和理论谱的比较。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

根据光子与物质的相互作用，用探测器测量单能光子，入射光子在进入探测器灵敏体积前可能会与屏蔽体、探测器外壳、空气等发生多种相互作用，进入探测器灵敏体积后与之发生的相互作用也不是单一的，而是包括多种不同的效应，因此单能光子的能量响应包括：全能峰、逃逸峰、反散射峰、康普顿坪、特征X射线等多种成分；而使用能谱仪测量X射线参考辐射得到的测量谱可以看作是能量连续的单能光子的能量响应的叠加，是脉冲高度谱，是连续谱X射线在谱仪中的沉积能量谱，为了得到连续谱X射线的注量谱，还需要进行能谱分析，常用的方法有：剥谱法、逆矩阵法、和反卷法等。

假设能谱仪对能量为E的单能光子的响应为R(E)，而辐射场的注量谱为F(E)，能谱仪测量连续谱得到的能谱为M(E)，则

$M\left(E\right)=\underset{0}{\overset{\∞}{\∫}}R(E,E_0)F\left(E_0\right)d{E}_0+\mathrm{\ϵ}---\left(1\right)$

上式中，ε表示由于探测器电子学系统所产生的随机噪声。将上式离散化得

$\left[\begin{array}{c}{M}_1\\ \·\\ \·\\ \·\\ M_m\end{array}\right]\left[\begin{array}{c}{R}_{11}{R}_{12}...R_{1n}\\ \begin{array}{cc}\·& \·\\ \·& \·\\ \·& \·\end{array}\\ R_{m1}{R}_{m2}...R_{\mathrm{mn}}\end{array}\right]\×\left[\begin{array}{c}{F}_1\\ \·\\ \·\\ \·\\ F_n\end{array}\right]+\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ \·\\ \·\\ \·\\ {\mathrm{\ϵ}}_2\end{array}\right]---\left(2\right)$

原则上上面的方程式，如果可以从探测器响应函数中分离出随机噪声项，可以用逆矩阵法来求解。但是一般很难将其分离出来，因此最终这部分响应被包含在探测器响应函数中，最终需求解的方程为

$\left[\begin{array}{c}{M}_1\\ \·\\ \·\\ \·\\ M_m\end{array}\right]\left[\begin{array}{c}{R}_{11}^{\′}{R}_{12}^{\′}...R_{1n}^{\′}\\ \begin{array}{cc}\·& \·\\ \·& \·\\ \·& \·\end{array}\\ R_{m1}^{\′}{R}_{m2}^{\′}...R_{\mathrm{mn}}^{\′}\end{array}\right]\×\left[\begin{array}{c}{F}_1\\ \·\\ \·\\ \·\\ F_n\end{array}\right]---\left(2\right)$

能谱测量和分析工作的工作流程可以分为两个部分：过滤X射线能谱的实验测量和能谱分析。能谱测量包括能谱测量系统建立，谱仪能量刻度，准直器设计。能谱分析包括能谱测量系统响应矩阵的确定和解谱程序的编写调试，响应矩阵的确定是采取蒙特卡罗模拟然后插值的方法，这就需要对探测器进行准确建模，要调整探测器模型中晶体的尺寸，要确定模拟中能谱的展宽特性等，最后达到模拟的单能光子在探测器中的能量沉积谱和实际测量得到的测量谱能很好符合。

本发明所提供的一种过滤X射线参考辐射的能谱分析方法，包括如下步骤：

(1)建立由探测器、数字化谱仪、计算机组成的能谱测量系统，并对能谱测量系统进行能量刻度，测量得到过滤X射线的测量谱；

(2)用蒙特卡罗模拟和数学插值的方法确定能谱测量系统的响应矩阵；

(3)对测量得到的测量谱进行剥谱法和反卷能谱法分析，得到了过滤X射线的注量谱，并与理论计算谱进行比较。

分析工作的基本工作流程如图1所示，下面结合具体工作流程对本发明进行详细的介绍。

●能谱测量系统的建立

在本发明中，选用的能谱测量系统包括：探测器、数字化谱仪、计算机。

探测器是选用适用于测量低能γ和X射线能谱的N型HPGe探测器，是由ORTEC公司生产的，型号为GMX25P4-70-S，已经在前面简单介绍，探测器的前窗为0.5mm铍窗加0.03mm的铝化聚脂膜，能量范围是3keV～2MeV，对1.33Mev的能量分辨率为1.9keV，工作电压为-4800V。探测器包括高纯锗晶体、前置放大器，电制冷装置。在能谱测量实验中，GB/T 12162中推荐的探测器离过滤X射线参考装置焦班距离为2.5m～3m，在实际测量中，选择在距焦斑3m处进行测量。

数字化谱仪是测量系统采集和处理的重要部件，选用的同样是ORTEC公司的产品，型号为DSPEC Plus^TM，它完成温度、高压状况的实时数据的采集、发出控制信号、报警，以实现HPGe探测器正常工作；完成死时间的计算与传输；完成脉冲信号的放大、成形；完成脉冲信号的采集(采用DSP芯片对信号进行实时采集)；完成电源的交直流转换、升压、稳压成比较稳定的高压电源以便输出到高压滤波模块中供探测器使用；完成高压直流向低压电源的转换，以便供系统中其他芯片使用；完成与计算机之间的数据传输与交流。

计算机实现状态数据的接受及处理、控制信号的输出、核数据的采集与处理、实时间的计算等一系列操作，使用的软件是数字化谱仪自带的软件MAESTRO32。

能谱测量系统的实验装置如图2所示，图中，1-X射线-射线管防护壳；2-初级光阑；3-快门；4-附加过滤片；5-限束光阑；6-准直器；7-铅屏蔽；8-高纯锗探测器。

能谱测量系统的能量刻度要求得到能量与道址的关系曲线，这条曲线代表谱曲线中的峰所对应的峰位道址与射线实际能量之间的关系。有了能量刻度曲线之后，寻峰程序在谱中找到的各个峰位就可以换算成所对应的能量。本发明采用的能谱测量系统中，解谱软件MAESTRO 32软件有设计好的能量刻度程序，应用该能量刻度程序，可以快速的进行能谱测量系统的能量刻度。

连接好能谱测量系统，加上探测器高压，调节上升时间，放大倍数等参数，确定了能谱测量系统工作时的各项参数之后，开始进行能量刻度。取一同位素刻度源，接近探测器放置，开始测量其能谱。该同位素刻度源发射的γ射线的能量应该已知，而且在能谱的测量范围之内。一段时间后，测量谱中出现明显的峰，将峰中心能量对应道址的能量改为同位素刻度源的能量。再另取一不同能量的同位素刻度源，用相同的办法进行能量刻度。反复进行几次，挑选的同位素刻度源能量要尽量覆盖整个测量范围，尤其是对即将测量连续X射线的能量范围，即300keV以下的能量范围。在本工作中采用的刻度源，有以下几种核素：¹³⁷Cs、²⁴¹Am、¹⁵²Eu、¹³³Ba。最后，通过验证，峰位道址所显示的能量和实际的入射光子的能量能很好的对应，结果列在表1中。

表1  能量刻度结果

在测量能谱时，如果谱仪的计数率太高，会造成探测器“堵塞”、死时间过高、出现和峰，为了减少入射到探测器晶体上的计数率、避免出现和峰，以及去除周围环境散射带来的复杂影响，设计了准直器(图2中的6所示)，起到屏蔽和准直的作用。准直器材料采用内部为铅，表面为2mm的不锈钢。在面向射束和垂直射束的方向铅屏蔽的厚度分别达到了7cm和5cm。5cm厚的铅，能将300keV的入射光子的强度衰减到10-10倍，所以该准直器的屏蔽效果已能都达到使用要求。准直器的准直孔开在“T”型的准直器上，孔径有φ2和φ5的，还有无开孔的实心堵头，用来测量本底。

根据准直器的设计分析，在测量能谱时，仍然会有部分反散射光子从探测器后部进入探测器。对此进行了实验验证，首先在关闭快门的情况下，测量本底的能谱，然后打开快门，把开孔的准直器换成不开孔的堵头，测量本底的能谱。结果显示，后者的能谱低能散射部分更多，与分析结果符合。因此在实际测量测量谱时，对每种参考辐射的测量谱扣除的本底谱，都应是在打开快门使用不开孔的堵头进行屏蔽后测量得到的。

●响应矩阵的建立

为了对测量谱进行剥谱法或反卷法能谱分析，必须要得到探测器在相应的测量条件下的响应矩阵(或全响应函数)。本发明中针对的过滤X射线的能量范围是6keV～300keV，所以要构建的响应矩阵是探测器对6keV～305eV的单能光子的响应矩阵。使用同位素源(¹³⁷Cs、¹³³Ba、²⁴¹Am等)，只能得到探测器对少数几个能量的探测效率及响应函数，而通过蒙特卡罗模拟计算的方法，可以得到多个探测器测量单能的光子探测效率及响应函数，这种方法简便、省时且准确度高。

在本发明的分析工作中得到响应函数所采用的的方法：用蒙特卡罗方法模拟探测器对同位素源光子的响应函数，在同样的实验条件下，用探测器实际测量同位素源得到实际的响应函数，两者相比较，调整蒙特卡罗方法模拟中探测器晶体的尺寸。用蒙特卡罗方法模拟探测器对多个单能光子的响应函数，以10keV为间隔，再用分段插值的方法得到其他能量的单能光子的响应，能量间隔1keV，构建响应矩阵。另外，在250keV～260keV之间的能量范围内，反散射峰和康普顿边缘相互交叠，用分段插值的方法会带来较大误差，因此，探测器对这各能量范围内的单能光子的响应函数，全部用蒙特卡罗模拟的办法得到。

蒙特卡罗方法又称随机抽样技巧或统计试验方法。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大的区别。它是以概率统计理论为基础的一种方法。蒙特卡罗的基本思想：当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。

由于蒙特卡罗方法在粒子输送方面的模拟中很有用，因此被广泛应用于实验核物理，反应堆物理，高能物理等领域。蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。

在本发明中，用蒙特卡罗方法解决的问题X光子在HPGe探测器中的能量沉积谱及响应函数的问题，首先对HPGe探测器及周围的其他实验装置和环境进行数学建模，然后通过模拟大量由放射源发射并通过HPGe探测器的单个X光子的历史，统计每个X光子在在HPGe探测器中沉积的能量，得到服从统计分布的物理量——能量在HPGe探测器中能量沉积的分布。

使用HPGe探测器测量单能光子得到的能谱，全能峰上计数随道数(即能量、脉冲振幅)的分布近似服从高斯分布，该高斯分布的特征参数依赖于探测器性能。为了使模拟计算得到的谱更接近实验测得的谱，得到更加准确的响应函数，所以模拟计算时也要考虑能谱的高斯展宽。

模拟计算中使用MCNP5程序提供的F8计数类型卡，记录光子在晶体中沉积的能量脉冲的幅度分布，然后选用计数特殊处理卡GEB对脉冲能量幅度分布进行高斯展宽。使用计数特殊处理卡GEB，需要设定三个参数，a、b和c，这三个参数设定后，可以计算出全能峰的半高宽，计算公式为：

$E_{\mathrm{FWHM}}=a+b\sqrt{E_r+c{E}_r^2}---\left(4\right)$

式中，E_r为入射γ射线的能量(单位为MeV)，E_FWHM为全能峰的半高宽(单位为MeV)，a参数的单位为MeV，b参数的单位为MeV^1/2，c参数单位为1。

首先，分别测量多个单能光子的能谱，分析得到的实验能谱，统计各个能量的单能光子对应的能量分辨率，即全能峰的半高宽，采用最小二乘法进行拟合，获得了参数a、b和c的值，带入公式计算各能量光子的能量分辨率与实测值进行比较。最后，确定当a＝0.0011344，b＝0.00003194，c＝610时，模拟计算得到的能量分辨率和实验得到的符合的很好，结果如下：

表2  实验和模拟得到的全能峰半高宽比较

为了确保蒙特卡罗方法模拟的探测器响应函数的精度，首先要确保对于单能光子的响应函数的正确性，这取决于对探测器及周围其他实验设备的准确建模。许多研究者发现使用蒙特卡罗方法计算HPGe探测器探测效率时，直接使用厂家提供的关于探测器组成结构的几何尺寸，特别是晶体的长度、半径和死层厚度等，得到的模拟计算结果与实际测量结果误差较大，进而发现厂家提供的探测器的结构几何尺寸与实际的尺寸有较大误差，对计算结果的准确性影响明显，为了使计算结果和测量结果相符合，需要精确确定探测器结构的几何尺寸。

确定探测器的结构几何尺寸，常用的方法有X射线成像、源扫描测量、实验解析等。在本工作中，通过实验测量和蒙特卡罗方法结合的方法确定HPGe探测器晶体的死层厚度、半径和长度。

首先，按照厂家提供的结构、尺寸对探测器内部晶体进行模拟。在不加准直器的情况下，把IAEA标准源放在不同方向确定位置的几个点上测量其能谱，并且用MCNP5模拟，比较模拟谱和测量谱的全能峰效率、峰康比等特征值，参考比较结果调整MCNP5模拟中探测器晶体的尺寸。

确定探测器晶体尺寸的实验如图3所示，图3中以O点为原点(0，0，0)，则g点坐标为(0，50，0)，i点坐标为(10，50，0)，i’点坐标为(-10，50，0)，a点坐标为(0，60，0)，d点坐标为(30，0，0)，单位为cm。

使用IAEA标准源¹³⁷Cs，调整晶体的长度和直径，因为¹³⁷Cs的能量(661.66keV)处于中能范围内，HPGe探测器对它的探测效率受死层厚度影响很小，并且在中能段可以看做是单能的，不会出现不同单能光子的响应互相叠加的情况；使用IAEA标准源²⁴¹Am调整晶体外部死层厚度，因为²⁴¹Am的能量(59.54keV)较低，HPGe探测器对它的探测效率受死层厚度影响很大。

厂家给出的探测器晶体长度为62.5mm，半径为26.5mm，外部死层厚度为0.3μm。按照这个晶体尺寸进行MCNP5的模拟计算，在图3中的a点、g点，对¹³⁷Cs的全能峰效率模拟计算结果比实际测量结果分别大13％和10％。

经过反复调节，最后确定HPGe探测器长度为58.5mm，半径为25.5mm，外部死层厚度为16.5μm，此时在各点对¹³7Cs和²⁴¹Am的全能峰效率模拟计算结果与实际测量结果相差小于2％。

确定能谱展宽特性以及完成探测器晶体尺寸调节后，模拟计算高纯锗X射线能谱仪对IAEA标准源¹³³Ba在没有屏蔽和准直器的情况下测量得到的能谱，并且与同样条件下实验测量得到的测量谱进行比较，测量谱归一到相对探测效率，结果如图4，图中纵坐标为归一后的计数率，横坐标为能量，单位为MeV。

从图4中可以发现对于¹³³Ba各个能量的全能峰，测量谱和模拟谱的结果能符合相当好，说明模拟谱的展宽特性、晶体尺寸调整的结果比较好。起初，在模拟能谱时没有考虑实验室的墙壁、地面、试验台等周围环境的其他物质，测量谱的康普顿坪明显高于模拟谱，结果分析，认为是周围环境引起的散射成分的影响，比如试验台、墙壁、地面等物质。之后，在模拟能谱时考虑了这些因素，对试验台、墙壁、地面等物质也都进行了简单的建模，模拟计算结果有了很大的改善，但是在对数坐标下依然能看出在康普顿边缘到全能峰之间的部分，测量谱的计数率略多于模拟谱的计数率，这可能是因为周围环境的模拟还是过于简单造成的。

对模拟得到的探测器响应函数，以如下方式进行分段内插，得到能量范围为6keV～305keV的光子在本工作中实际使用的能谱测量系统中的完整的响应矩阵。

假设有两个探测器响应函数R₁(E₁，e)，R₂(E₂，e)，根据能谱中的特征能点，在能量尺度上，将能谱分为4个区域：(5，E_b)，(E_b，E_C)，(E_C，E_C-10)，(E_C-10，E_P+10)。其中，E_b＝E_P/(1+2E_P/m_ec²)，E_C＝E_P-E_b及E_P分别为反散射峰，康普顿边缘及光电峰所对应的能量。对于内插得到的响应函数R_i(E_i，e)，

首先确定对应于R_i(E_i，e)的分割点

$E_b^{\′}=\frac{E_i-E_2}{E_1-E_2}[E_b^1-E_b^2]+E_b^2---\left(5\right)$

$E_C^{\′}=\frac{E_i-E_2}{E_1-E_2}[E_C^1-E_C^2]+E_C^2---\left(6\right)$

$E_P^{\′}=\frac{E_i-E_2}{E_1-E_2}[E_P^1-E_P^2]+E_P^2---\left(7\right)$

对于R_i(E_i，e)每一个分区中的每一个能量点，分别在R₁(E₁，e)，R₂(E₂，e)中确定与之相对应的能量点。例如对于e_i∈(E′_b，E′_C)

$e_{j1}=\frac{e_j-E_b^{\′}}{E_C^{\′}-E_b^{\′}}[E_C^1-E_b^1]+E_b^1---\left(8\right)$

$e_{j2}=\frac{e_j-E_b^{\′}}{E_C^{\′}-E_b^{\′}}[E_C^2-E_b^2]+E_b^2---\left(9\right)$

在响应函数R₁(E₁，e)和R₂(E₂，e)中，确定分别对应于e₁和e₂的响应值，并由此确定R_i(E_i，e)中e_i所对应的响应值

$R_i(E_i,e_j)=\frac{E_i-E_2}{E_1-E_2}[R_1(E_1,e_{j1})-R_2(E_2,e_{j2})]+R_2(E_2,e_{j2})---\left(10\right)$

由上述分段内插算法最终得到的响应矩阵。

●从测量谱得到注量谱的分析

从测量谱得到注量谱，就是解式(3)中方程的问题。众多的研究结果表明，由于式(3)中的方程，包含了随机噪声，并且测量谱及响应函数计算本身亦包含了统计涨落，因此直接解出来的结果存在着较大的震荡，因此目前一般都采用迭代法求解。

反卷过程，从数学看就是解一个线性方程组M＝R×F，F为注量谱，R为响应矩阵，M是实验测量得到的测量谱；解这个线性方程组的问题是一个典型的病态问题；对于X射线参考辐射场能谱的测量及根据探测器响应函数的反卷，已有不少可以借鉴的算法。在这里采用一种性能较好的反卷算法：Gold算法。这个算法其实是来源于Van Cittert迭代算法(实际上其可看成是多元函数极值搜索算法中的梯度搜索方法)。只不过，在每一步迭代过程中，其采用了一个可变的松弛因子。设测量谱向量为M，探测系统响应矩阵为R_M×N，入射辐射注量谱为F，则所要求解的问题是一个线性方程组M＝R·F。若再对其进行如下的操作。

对于上面最终得到的方程，Van Cittert迭代算法为

F⁽ⁿ⁺¹⁾＝F⁽ⁿ⁾+μ(Z-AF⁽ⁿ⁾)    (12)

上式中，n是迭代次数，μ可看成是梯度搜索算法中的搜索步长。如果将μ选为一个局部可变的松弛因子

${\mathrm{\μ}}_i=\frac{F_j^{\left(n\right)}}{{\mathrm{\Σ}}_{j=1}^N{A}_{\mathrm{ij}}{F}_j^{\left(n\right)}}---\left(13\right)$

并将其带入上式中，则得到最终的Gold算法

$F^{(n+1)}=\frac{Z_i{F}_i^{\left(n\right)}}{d_i},i=<1,n>---\left(14\right)$

上式中，d＝A·F⁽ⁿ⁾。

使用MATLAB7编写解谱程序，MATLAB既是一种直观、高效的计算机语言，同时又是一个科学计算平台。解谱程序流程如图5所示。“能道修正”部分的主要作用是重新划分能道，并把原能道的计数合理的分配到新的能道中。测量谱中310.09keV对应的能道为1297，每个能道的能量间隔大约为0.24keV。重新划分能道以后，每个能道的能量间隔为1keV。经过“能道修正”后，重建的测量谱，包括300道，能量范围：5.5keV～305.5keV，各道中心能量为(6，7，......，305)。

在本发明中用Gold算法分析得到的各参考辐射的注量谱与理论计算得到的各参考辐射的注量谱，按连续谱峰值归一后一一进行比较。图6和图7分别为低空气比释动能率系列和窄谱系列，解谱得到的注量谱和理论谱的比较(都是经过归一的)。

从图6和图7中注量谱和理论谱的比较可以看出，窄谱系列300kV、250kV、200kV和低空气比释动能率系列240kV、210kV五个参考辐射的注量谱和理论计算谱符合的不好，主要表现为在注量谱的峰左边有一个拖尾(或者说平台)，这可能是因为一部分入射X射线在准直器内发生康普顿散射造成的。

窄谱系列20kV、15kV两种参考辐射的注量谱和理论计算谱符合的也不是很好，这是因为目前分析注量谱和理论计算谱采用的能道宽度都为1keV，相对于最大能量为20keV或15keV的这些能量较低的参考辐射来说，能道宽度太宽，无法很好的表征出参考辐射的谱特性。

得到了过滤X射线的注量谱，计算过滤X射线参考辐射的一些特征参数，并与GB/T 12162给出的推荐值进行比较，以此来从另一个角度验证本发明得到的注量谱的准确性。

平均能量

对于连续谱X射线，平均能量的定义是由下列公式定义的：

$\stackrel{\&OverBar;}{E}=\frac{{\&Integral;}_{E_{\min }}^{E_{\max }}E{\mathrm{\&Phi;}}_E\left(E\right)\mathrm{dE}}{{\&Integral;}_{E_{\min }}^{E_{\max }}{\mathrm{\&Phi;}}_E\left(E\right)\mathrm{dE}}---\left(15\right)$

式中：Φ_E是Φ(E)能量对E的微商，Φ(E)是能量在O至E之间的光子总注量。

实际计算时，得到的注量谱为离散的，所以采用下面的公式：

$\stackrel{\&OverBar;}{E}=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i{\mathrm{\&Phi;}}_{E_i}}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&Phi;}}_{E_i}}---\left(16\right)$

能量分辨率

谱的能量分辨率(半峰高处的全高宽)是由下列公式定义的比值，以百分数表示：

$R_E=\frac{\mathrm{\&Delta;E}}{E}\&times;100\%---\left(17\right)$

式中：ΔE——对应于谱的最大纵坐标一半处的谱宽度。注：在谱中存在荧光辐射的情况下，测量的谱宽度仅仅根据连续谱求出。

空气比释动能的半值层

半值层(空气比释动能)HVL或HVLX的定义，是使辐射束的空气比释动能率减少到其初始值的一半的指定材料的厚度。确定半值层可以用实验的方法或者计算的方法。

同质系数

同质系数是指第一半值层(1stHVL)与第二半值层(2stHVL)(空气比释动能)的比值：

$h=\frac{1\mathrm{stHVL}}{2\mathrm{ndHVL}}---\left(18\right)$

表3  平均能量计算值和GB/T 12162推荐值的比较

注：表中参考辐射的表示方法为字母加数字，字母代表过滤X射线的系列，数字代表X光机的管电压，如“N15”指的是15kV窄谱系列，“L240”指的是240kV低空气比释动能率系列。

1)半值层材料为铝(Al)；

2)半值层材料为铜(Cu)。

剂量转换系数

在GB/T 12162中给出了对于各种参考辐射(不同能量的单能光子，不同系列的连续谱过滤X射线)，从空气比释动能到周围剂量当量、个人剂量当量的转换系数。在得到连续谱过滤X射线的注量谱后，可以利用公式(19)计算从空气比释动能到周围剂量当量、个人剂量当量的转换系数。

$h_{\mathrm{pK}}=\frac{{\&Integral;}_{E_{\min }}^{E_{\max }}\frac{d{K}_a\left(E\right)}{\mathrm{dE}}{h}_{\mathrm{pK}}\mathrm{dE}}{{\&Integral;}_{E_{\min }}^{E_{\max }}\frac{d{K}_a\left(E\right)}{\mathrm{dE}}\mathrm{dE}}---\left(19\right)$

式中h_pK是对于单能光子的转换系数；dK_a(E)/dK即每一能道对应的比释动能。

表4  剂量转换系数的计算值、GB/T 12162推荐值的比较

注：表中参考辐射的表示方法为字母加数字，字母代表过滤X射线的系列，数字代表X光机的管电压，如“N15”指的是15kV窄谱系列，“L240”指的是240kV低空气比释动能率系列。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (8)

1.一种过滤X射线参考辐射的能谱分析方法，包括如下步骤：

(1)建立由探测器、数字化谱仪、计算机组成的能谱测量系统，并对能谱测量系统进行能量刻度，测量得到过滤X射线的测量谱；

(2)用蒙特卡罗模拟和数学插值的方法确定能谱测量系统的响应矩阵；

(3)对测量得到的测量谱进行剥谱法和反卷能谱法分析，得到了过滤X射线的注量谱，并与理论计算谱进行比较。

2.如权利要求1所述的过滤X射线参考辐射的能谱分析方法，其特征在于：所述的能谱测量系统的探测器为HPGe探测器。

3.如权利要求1或2所述的过滤X射线参考辐射的能谱分析方法，其特征在于：步骤(1)中对能谱测量系统进行能量刻度的方法如下：取一同位素刻度源，接近探测器放置，开始测量其能谱；一段时间后，测量谱中出现明显的峰，将峰中心能量对应道址的能量改为该同位素刻度源的能量；再另取一不同能量的同位素刻度源，用相同的办法进行能量刻度；反复进行几次，挑选的同位素刻度源能量应能够覆盖即将测量连续X射线的能量范围。

4.如权利要求3所述的过滤X射线参考辐射的能谱分析方法，其特征在于：步骤(1)中对能谱测量系统进行能量刻度所选用的同位素刻度源包括：¹³⁷Cs、²⁴¹Am、¹⁵²Eu、¹³³Ba。

5.如权利要求1或2所述的过滤X射线参考辐射的能谱分析方法，其特征在于：步骤(2)中所述的蒙特卡罗模拟方法是首先对HPGe探测器及周围的其他实验装置和环境进行数学建模，然后通过模拟大量由放射源发射并通过HPGe探测器的单个X光子的历史，统计每个X光子在HPGe探测器中沉积的能量，得到服从统计分布的在HPGe探测器中能量沉积的分布。

6.如权利要求5所述的过滤X射线参考辐射的能谱分析方法，其特征在于：在蒙特卡罗模拟计算中使用MCNP5程序提供的F8计数类型卡记录光子在晶体中沉积的能量脉冲的幅度分布，然后选用计数特殊处理卡GEB对脉冲能量幅度分布进行高斯展宽。

7.如权利要求5所述的过滤X射线参考辐射的能谱分析方法，其特征在于：步骤(2)中所述的数学插值的方法是采用分段内插的方式，得到能量范围为6keV～305keV的光子在能谱测量系统中的完整的响应矩阵。

8.如权利要求1或2所述的过滤X射线参考辐射的能谱分析方法，其特征在于：步骤(3)中从测量谱得到注量谱是通过Gold反卷算法求解线性方程组M＝R×F，其中，F为注量谱，R为响应矩阵，M是实验测量得到的测量谱。

Images