CN105785424A - 一种碲锌镉探测器伽玛谱全能峰非线性拟合算法 - Google Patents

Abstract

本发明涉及一种碲锌镉探测器伽玛谱全能峰非线性拟合算法。该算法提供了碲锌镉探测器伽玛能谱全能峰拟合函数，在拟合函数的非线性拟合过程中，借鉴工业上物理退火的物理过程和Metropolis准则，结合最小二乘Levenberg‑Marquardt(LM)迭代算法，采用内、外两个循环，进行非线性拟合，既提高了迭代计算的速度，又达到全局收敛的目的，克服了迭代结果对拟合参数初始值过度依赖的缺点。

Description

一种碲锌镉探测器伽玛谱全能峰非线性拟合算法

技术领域

本发明涉及一种针对碲锌镉探测器伽玛能谱全能峰非线性拟合的计算方法，属于核辐射监测技术领域。

背景技术

碲锌镉(CZT)探测器属于化合物半导体探测器，可在常温下使用，体积小，探测装置简单，适于便携式测量。与其它半导体探测器以及碘化钠(NaI)探测器相比，CZT探测器的显著优点是具有宽的禁带和低的电离能，这使它在常温下具有很好的能量分辨率，同时高的原子序数提高了全能峰的本征效率。随着CZT晶体制备技术的提高以及探测器电极的特殊设计方法的使用，晶体的尺寸及其性能不断得到提高和改善。目前用于科学研究的CZT探测器能量分辨率达到1％(662keV)，商用探测器能量分辨率达到2％左右，明显优于NaI探测器。自从CZT探测器出现以来，在国家安防、空间物理、医学成像以及工业等领域迅速得到应用。

CZT探测器测得的伽玛能谱的全能峰最显著的特点是存在低能尾巴，在高能段尤为严重。这是因为CZT晶体中载流子平均漂移长度偏低，电荷收集过程中空穴很容易被陷阱俘获而不能到达电极造成的，同时导致能量分辨率下降。如果直接套用高纯锗(HPGe)或者NaI探测器伽玛谱的解谱算法会产生较大误差，因此需要根据CZT探测器伽玛谱自身的特点研究与其相适应的拟合函数。

1999年出版的《NuclearInstrumentandMethodinPhysics》第A422卷期刊中《Cadmiumzinctelluridespectralmodeling》(159～163页)一文中给出了CZT探测器测得的伽玛谱全能峰拟合函数。2006年清华大学工程物理系艾宪芸的博士论文《碲锌镉探测器性能分析及其γ谱解析方法研究》中进一步验证了这种拟合函数的适用性。这种拟合函数能较好的描述碲锌镉探测器所测的伽玛谱光电全能峰。但这种拟合函数中共有10个未知参数，需要对其进行非线性拟合计算。

在伽玛能谱的非线性拟合计算中通常采用最小二乘拟合方法，最小二乘拟合方法的关键是寻找最佳解使其与真值的均方差最小，最常用的计算方法是Levenberg-Marquardt(LM)法，LM算法是在高斯-牛顿法的基础上融合了梯度下降法后的改进方法。目前已经成为非线性最小二乘拟合算法的标准算法。

使用LM算法对伽玛谱全能峰进行非线性拟合，拟合效果较好，但存在如下缺陷：拟合结果对非线性参数初始值的依赖性较强，有时甚至得不到正确结果，原因是LM算法不是全局收敛的。

发明内容

本发明的目的针对上述存在的问题，提供了一种LM算法的改进算法，借鉴工业上模拟退火的物理过程和Metropolis准则，在LM算法中，改变迭代过程中的搜索方向，即不仅可以向目标函数增大的方向搜索，也能向目标函数减小的方向搜索，从而可以从局部极值中爬出，不会陷在局部极值中，达到全局收敛的目的，同时也克服了迭代结果对拟合参数初始值过度依赖的缺点。提供一种碲锌镉探测器伽玛谱全能峰非线性拟合算法。该方法也适用于其他探测器所测伽玛谱全能峰的非线性拟合计算。

本发明的碲锌镉探测器伽玛谱全能峰非线性拟合算法具体步骤如下：

第一步：确定CZT探测器伽玛谱全能峰的拟合函数式为：

f(i，p)＝G(i)+B(i)+S(i)+D(i)

式中，G(i)＝H_g×exp[-(i-i₀)²/(2σ²)]；

B(i)＝A₅+A₆×i；

S(i)＝H_s×H_g×erfc[(i-i₀)/(σ2^1/2)]；

$\begin{array}{c}D\left(i\right)=H_t\×H_g\×\exp \[(i-i_0)/\left(T_s\mathrm{\σ}\right)\]\×erfc\[(i-i_0)/\left(\mathrm{\σ}{2}^{1/2}\right)\\ \begin{array}{cc}+1/\left(T_s{2}^{1/2}\right)\];& erfc\left(x\right)=2/\sqrt{\mathrm{\π}}{\∫}_x^{\mathrm{\∞}}\exp (-t^2)dt\end{array}\end{array};$

式中：H_s为本底台阶高度；H_g为高斯峰的峰高；i_o＝A₁+A₂E_peak为峰位；σ²＝A₃+A₄E_peak为方差；H_t为尾巴高度；T_s为指数尾巴的反斜率；i是道址；

第二步：确定需要拟合的参数共10个：H_s、H_g、A₁、A₂、A₃、A₄、A₅、A₆、H_t、T_s；

第三步：设定CZT探测器伽玛谱全能峰的测量谱为y(i)，i表示道址；

第四步：设定每个拟合参数的取值范围，其极大、极小值分别存入数组AH和AL中；

第五步：温度初始化，对每个拟合参数在其取值范围内均匀抽样得到一组拟合参数的初始值，该初始参数向量记作p₁，把p₁带入拟合函数得到f(i，p₁)，计算残差平方和：重复该过程n次，其中n大于10，得到一组χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²值，选择其中最大的数值作为初始温度T₀＝max(χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²)；

第六步：外循环计算，拟合参数向量为p_i，对应的温度为T_i＝χ_i ²，目标函数值为进行第七步的内循环计算，得到一组新的参数向量p_j，计算目标函数值计算Δ_E＝E(T_j)-E(T_i)，如果Δ_E＜0，则接受T_j为新状态，温度更新为T_j＝χ_j ²，否则，计算几率p＝exp(-Δ_E/T_i)，在[0，1]区间上均匀抽样，得到抽样值c，如果c小于p，则接受T_j为新状态，温度更新为T_j＝χ_j ²，否则，保持T_i状态不变，重复第七步的内循环计算；

第七步：内循环计算，输入量为拟合参数向量p_i，拟合函数是f(i，p_i)，计算其雅可比矩阵计算拟合函数与测量谱数据的差值：ε＝y(i)-f(i，p_i)，新的拟合参数向量取值为：p_j＝p_i+δ_p，δ_p的取值由下式计算获得：J^TJδ_p＝J^Tε；

第八步：当达到理想的低温时，即得到最优解时，外循环计算结束；外循环计算结束判定条件：|Δ_E|小于事先设定的极小值。

初始温度值采用一组残差平方和数据中的最大值；迭代过程中，当温度或者残差平方和降低，则把减小后的残差平方和作为新一轮迭代的温度值；当温度升高，则根据Metropolis准则，判断是否接受升温后的状态为新状态，不接受则保持当前温度状态不变。

拟合算法采用模拟退火和LM相结合的方法，外循环采用模拟退火方法，内循环采用LM方法。

工业上物理退火的一般过程是：设粒子在温度T时的状态为i，对应的系统能量是E_i，然后假设粒子在一个随机扰动作用下进入一个新的状态j，此时对应的系统能量是E_j，如果E_j＜E_i，则接受新状态j为当前状态，否则，根据粒子处于新状态j的几率来判断是否接受新状态j为当前状态，如不接受新状态，则保持原来的i状态，然后进行上述循环迭代，直到目标函数满足预期要求。

对于碲锌镉探测器伽玛谱全能峰非线性拟合过程，初始温度T₀的设定如下：对所有需要拟合的未知参数随机抽样得到一组参数值，把这组参数值代入拟合函数，然后求出拟合函数与伽玛实测谱的残差平方和χ₁ ²，重复该过程n次，得到一组χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²值，选择其中最大的数值作为初始温度T_i＝T₀。

在温度T_i状态下，目标函数值为E(T_i)，利用LM算法得到新的拟合参数，计算拟合函数与伽玛实测谱的残差平方和χ_i+1 ²，此时的目标函数值为E(T_i+1)，计算Δ_E＝E(T_i+1)-E(T_i)，如果Δ_E＜0，则接受T_i+1为新状态，温度T_i+1＝χ_i+1 ²，否则，计算几率：p＝exp(-Δ_E/T(i))，在[0，1]区间上均匀抽样，得到c，如果c小于p，则接受T_i+1为新状态，温度T_i+1＝χ_i+1 ²，否则，保持T_i状态不变，重复迭代计算。

在上述迭代过程中，有两个循环，外循环由系统温度T_i控制，其作用是在总体保证目标函数值下降的过程中，允许目标函数值在一定范围内上升跳变，从而使迭代过程能跳出局部极小值，达到全局极小值的目的。内循环采用LM算法，LM算法的特点是非等步长迭代，即距离真值较远时，迭代步长取小，虽然收敛速度较慢，但能保证局部收敛，当接近真值时，迭代步长增加，加速收敛速度。

本发明的有益效果：结合了LM算法和模拟退火算法的优点，既提高了迭代计算的速度，又达到全局收敛的目的，克服了迭代结果对拟合参数初始值过度依赖的缺点。

具体实施方式

下面结合附图和实施方式对本发明的发明内容作进一步说明。

第一步：确定CZT探测器伽玛谱全能峰的拟合函数如(1)式：

f(i，p)＝G(i)+B(i)+S(i)+D(i)(1)

其中，G(i)＝H_g×exp[-(i-i₀)²/(2σ²)]；

B(i)＝A₅+A₆×i；

S(i)＝H_s×H_g×erfc[(i-i₀)/(σ2^1/2)]；

$\begin{array}{c}D\left(i\right)=H_t\×H_g\×\exp \[(i-i_0)/\left(T_s\mathrm{\σ}\right)\]\×erfc\[(i-i_0)/\left(\mathrm{\σ}{2}^{1/2}\right)\\ \begin{array}{cc}+1/\left(T_s{2}^{1/2}\right)\];& erfc\left(x\right)=2/\sqrt{\mathrm{\π}}{\∫}_x^{\mathrm{\∞}}\exp (-t^2)dt\end{array}\end{array};$

其中，H_s为本底台阶高度；

H_g为高斯峰的峰高；

i_o＝A₁+A₂E_peak为峰位；

σ²＝A₃+A₄E_peak为方差；

H_t为尾巴高度；

T_s为指数尾巴的反斜率；

i是道址。

第二步：由(1)式知，需要拟合的参数共10个：H_s、H_g、A₁、A₂、A₃、A₄、A₅、A₆、H_t、T_s，其中H_s、H_g、A₅、A₆、H_t为线性拟合参数，A₁、A₂、A₃、A₄、T_s为非线性拟合参数。

第三步：设定CZT探测器伽玛谱全能峰的测量谱为y(i)，i为道址；

第四步：设定每个拟合参数的取值范围，其极大、极小值分别存入数组AH和AL中，即AH＝[H_s-min、H_g-min、A_1-min、A_2-min、A_3-min、A_4-min、A_5-min、A_6-min、H_t-min、T_s-min]，AL＝[H_s-max、H_g-max、A_1-max、A_2-max、A_3-max、A_4-max、A_5-max、A_6-max、H_t-max、T_s-max]；

第五步：温度初始化：对每个拟合参数在其取值范围内均匀抽样得到一组拟合参数的初始值，该初始参数向量记作p₁，把p₁带入拟合函数得到f(i，p₁)，计算残差平方和：重复该过程n次(n＞10)，得到一组χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²值，选择其中最大的数值作为初始温度：T₀＝max(χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²)。

第六步：外循环计算：拟合参数向量为p_i，对应的温度为T_i＝χ_i ²，目标函数值为进行S7内循环计算，得到一组新的参数向量p_j，计算目标函数值为计算Δ_E＝E(T_j)-E(T_i)，如果Δ_E＜0，则接受T_j为新状态，温度更新为T_j＝χ_j ²，否则，计算几率：p＝exp(-Δ_E/T_i)，在[0，1]区间上均匀抽样，得到抽样值c，如果c小于p，则接受T_j为新状态，温度更新为T_j＝χ_j ²，否则，保持T_i状态不变，重复S7内循环计算。

第七步：内循环计算：内循环开始的输入量为拟合参数向量p_i，拟合函数是f(i，p_i)，计算其雅可比矩阵：拟合函数与测量谱数据的差值ε＝y(i)-f(i，p_i)，新的拟合参数向量取值为：p_j＝p_i+δ_p，δ_p的取值由下式计算获得：J^TJδ_p＝J^Tε。

第八步：当达到最低温度，即得到最优解后，外循环计算结束，外循环计算结束判定条件：|Δ_E|小于事先设定的某极小值。

碲锌镉探测器伽玛谱全能峰非线性拟合算法。利用工业上物理退火的原理，结合最小二乘Levenberg-Marquardt(LM)迭代算法，采用内、外两个循环，对碲锌镉探测器伽玛谱全能峰的拟合函数进行了非线性拟合，较好得解决了多个非线性参数的拟合问题，该方法既提高了迭代计算的速度，又达到全局收敛的目的，克服了迭代结果对拟合参数初始值过度依赖的缺点。

Claims (1)

1.一种碲锌镉CZT探测器伽玛全能峰非线性拟合算法，其特征在于该拟合算法具体步骤如下：

第一步：确定CZT探测器伽玛谱全能峰的拟合函数式为：

f(i，p)＝G(i)+B(i)+S(i)+D(i)

式中：G(i)＝H_g×exp[-(i-i₀)²/(2σ²)]；

B(i)＝A₅+A₆×i；

S(i)＝H_s×H_g×erfc[(i-i₀)/(σ2^1/2)]；

D(i)＝H_t×H_g×exp[(i-i₀)/(T_sσ)]×erfc[(i-i₀)/(σ2^1/2)+1/(T_s2^1/2)]；

$\mathrm{erfc}\left(x\right)=2/\sqrt{\mathrm{\π}}{\∫}_x^{\∞}\exp \left({-t}^2\right)\mathrm{dt};$

式中：H_s为本底台阶高度；H_g为高斯峰的峰高；i_o＝A₁+A₂E_peak为峰位；σ²＝A₃+A₄E_peak为方差；H_t为尾巴高度；T_s为指数尾巴的反斜率；i为道址；

第二步：确定需要拟合的参数共10个：H_s、H_g、A₁、A₂、A₃、A₄、A₅、A₆、H_t、T_s；

第三步：设定CZT探测器伽玛谱全能峰的测量谱为y(i)，i表示道址；

第四步：设定每个拟合参数的取值范围，其极大、极小值分别存入数组AH和AL中；

第五步：温度初始化，对每个拟合参数在其取值范围内均匀抽样得到一组拟合参数的初始值，该初始参数向量记作p₁，把p₁带入拟合函数得到f(i，p₁)，计算残差平方和：重复该过程n次，其中n大于10，得到一组χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²值，选择其中最大的数值作为初始温度T₀＝max(χ₁ ²、χ₂ ²、χ₃ ²、...、χ_n ²)；

第六步：外循环计算，拟合参数向量为p_i，对应的温度为T_i＝χ_i ²，目标函数值为进行第七步的内循环计算，得到一组新的参数向量p_j，计算目标函数值计算Δ_E＝E(T_j)-E(T_i)，如果Δ_E＜0，则接受T_j为新状态，温度更新为T_j＝χ_j ²，否则，计算几率p＝exp(-Δ_E/T_i)，在[0，1]区间上均匀抽样，得到抽样值c，如果c小于p，则接受T_j为新状态，温度更新为T_j＝χ_j ²，否则，保持Ti状态不变，重复第七步的内循环计算；

第七步：内循环计算，输入量为拟合参数向量p_i，拟合函数是f(i，p_i)，计算其雅可比矩阵计算拟合函数与测量谱数据的差值：ε＝y(i)-f(i，p_i)，新的拟合参数向量取值为：p_j＝p_i+δ_p，δ_p的取值由下式计算获得：J^TJδ_p＝J^Tε；

第八步：当达到理想的低温时，即得到最优解时，外循环计算结束；外循环计算结束判定条件：|Δ_E|小于事先设定的极小值；

初始温度值采用一组残差平方和数据中的最大值；迭代过程中，当温度或者残差平方和降低，则把减小后的残差平方和作为新一轮迭代的温度值；当温度升高，则根据Metropolis准则，判断是否接受升温后的状态为新状态，不接受则保持当前温度状态不变；

拟合算法采用模拟退火和LM相结合的方法、外循环采用模拟退火方法、内循环采用LM方法。