UP | HOME

SPAN5 中的 γ 谱的解谱技术

目录

SPAN 是原子能量的 WANG Liyu 针对 DOS 开发的一款 γ 解谱软件。广泛应用在原子能院内部的相关工作当中。目前仍然在工作中得到较广泛的应用。此处摘录该软件手册中介绍的 γ 解谱办法,仅供参考。

这部分内容来自于 "SPAN V5.1" 软件手册的第 15 章。

活度计算

库中的数据与谱中的某个峰匹配的话,则该核素特定的活度计算如下:

\begin{equation*} A = \frac{N_{t}}{\gamma \varepsilon T_{c}} \cdot \frac{\lambda T_{t}}{1 - e^{-\lambda T_{t}}} \end{equation*}

此处:

  • \(A\) – 活度 (Bq)
  • \(N_{t}\) – 净峰面积(所考虑的光电峰计数)
  • \(\lambda\) – 衰变常数 (\(s^{-1}\))
  • \(T_{c}\) – 采集活时间 (s)
  • \(T_{t}\) – 采集实时间 (s)
  • γ – 感兴趣 γ 射线的 γ 分支比 (每次退激中出射特定能量的 γ 射线的数量)
  • \(\varepsilon\) – 感兴趣能区的探测效率

活度增长和衰变修正

活化的靶核放射性活度会增长,则放射性同位素的活度为:

\begin{equation*} A = N \Phi \sigma \left( 1- e^{-\lambda T_{i}} \right) e^{-\lambda T_{d}} \end{equation*}

此处:

  • \(A\) – 活度 (Bq)
  • \(N\) – 靶核数量
  • \(\Phi\) – 中子通量 (\(\mbox{s}^{-1}\cdot \mbox{cm}^{-2}\))
  • \(\sigma\) – 中子活化截面 ( cm^2 )
  • \(\lambda\) – 衰变常数 (s-1)
  • \(T_{i}\) – 活化时间 (s)
  • \(T_{d}\) – 衰变时间 (s)
  • \(N = 6.023 \times 10^{23} \times m \times C \times \rho / A_t\)
    • \(m\) – 样品质量 (g)
    • \(C\) – 元素浓缩度
    • \(\rho\) – 靶同伴素在元素中的自然含量丰度
    • \(A_{t}\) – 原子质量 (g)

本底

使用积分法峰处理技术,可生成一个阶段的本底曲线,其形状趋近峰下面的本底,从而可以从谱中被直接扣除。该技术计算可接受的净峰面积。其原理如后。

γ 谱中峰下的本底是阶梯状的,在很小区间内,其步长与峰高是成比例的,因此可以假设第 i 道的本底步长与该道的净计数成正比:

\begin{equation*} dH = \frac{\left(A_{i}-B_{i}\right) \cdot H_{t}}{S} \end{equation*}

此处:

  • \(dH\) – 第 i 道的本底步长
  • \(A_i\) – 第 i 道计数
  • \(B_i\) – 第 i 道本底
  • \(S\) – 净峰面积, \(S = \sum\limits^{r}_{i=l}\left(A_i-B_i\right)\)
    • \(l\) – 峰的左边界
    • \(r\) – 峰的右边界

总的本底步长, \(H_t\), 是峰前后的本底水平的差值(参见如后方程)。

\begin{eqnarray*} H_t & = & a_1 + b_1 \cdot c - a_2 - b_2 \cdot c \\ c & = & \frac{\sum\limits^{r}_{i=l}\left[ \left(A_i - B_i\right) \cdot i \right]}{\sum\limits^r_{i=l}\left(A_i - B_i\right)} \\ \end{eqnarray*}

这些本底可以定义为两条直线,如下:

\begin{eqnarray*} y_1 &= a_1 + b_1 x &\quad\quad \mbox{left side} \\ y_2 &= a_2 + b_2 x &\quad\quad \mbox{right side} \\ \end{eqnarray*}

则积分后的峰本底 \(B_i\) 定义如下:

\begin{equation*} B_i = a_1 + b_1 \cdot i - \frac{H_t}{S} \times \left\{ \left[ \sum\limits^i_{k=l}\left(A_k - B_k\right) \right] - \frac{A_i - B_i}{2} \right\} \end{equation*}

当然, \(dH\), \(B_i\), \(c\) 和 \(H_t\) 是相互关联的,但并非绝对。\(B_i\) 的初值由前面的两条直线给出,最终的 \(B_i\) 由渐近法得到。

符合叠加修正

如果在 γ 射线谱仪的分辨时间内从放射性核依次连续发射两个或多个光子,则会发生符合叠加。比如, γ_1 和 γ_2 的能量分别为 E_1 和 E_2 , 从谱仪分辨时间内发射。没有符合叠加时, γ_1 的全能峰脉冲发射率为:

\begin{equation*} n_{10} = A \times \rho_1 \times \varepsilon_1 \end{equation*}

此处:

  • \(A\) – 源活度
  • \(\rho_1\) – γ_1 的发射几率
  • \(\varepsilon_1\) – 在能量 E_1 处的全能峰效率

实际上,实验观测到的全能峰脉冲的发射率为 \(n_1\), 它将比 \(n_{10}\) 要小。由于每条 γ_1 紧跟一条符合的 γ_2 , 则导致两条 γ 射线都被探测到,从而叠加成一个脉冲,其能量与任何一条原始脉冲的值都不相同。探测 γ_2 的几率与总的探测效率 \(\varepsilon_t\) 相等,则有:

\begin{equation*} n_1 = A \cdot \rho_1 \cdot \varepsilon_1 - A \cdot \rho_1 \cdot \varepsilon_1 \cdot \varepsilon_{t2} = A \cdot \rho_1 \cdot \varepsilon_1 \cdot \left( 1 - \varepsilon_{t2} \right) \end{equation*}

则 γ_1 的修正系数为:

\begin{equation*} C_1 = \frac{n_{10}}{n_1} = \frac{1}{1-\varepsilon_{t2}} \end{equation*}

亦可得到 γ_2 的修正系数:

\begin{equation*} C_2 = \frac{1}{1 - \rho_1 \cdot \varepsilon_{t1}/\rho_2} \end{equation*}

此处:

  • \(\rho_2\) – γ_2 的出射几率
  • \(\varepsilon_{t1}\) – γ_1 的总效率

有些核的 γ 衰变链和直接关联很复杂,则可使用因子 \(w_i\) 来修正 \(\varepsilon_{t1}\) 和 \(\varepsilon_{t2}\), 则 \(C_1\) 和 \(C_2\) 表达式为:

\begin{eqnarray*} C_1 & = & \frac{1}{1 - w_2 \cdot \varepsilon_{t2}} \\ C_2 & = & \frac{1}{1-w_1 \cdot \rho_1 \cdot \varepsilon_{t1} / \rho_2} \\ \end{eqnarray*}

SPAN 软件中每条衰变链中只允许标识两条 γ 射线,当 \(\rho_1 \ge \rho_2\) 时使用上面方程组中的第一个方程,反之则使用第二个。而因子 \(w_1\) 和 \(w_2\) 由实验获得。

效率刻度曲线

探测器在特定能量的探测效率可由下式表示:

\begin{equation*} \log(\eta) = \begin{cases} a_1 + b_1 \log(E) + c_1 \log^2(E) & E > E_{m1} \\ a_2 + b_2 \log(E) + c_2 \log^2(E) & E \le E_{m1} \quad\text{or}\quad E>E_{m2} \\ a_3 + b_3 \log(E) + c_3 \log^2(E) & E \le E_{m2} \end{cases} \end{equation*}

此处:

  • \(\eta\) – 探测效率
  • \(E\) – 能量
  • \(E_{mi}\) – 中点能量,并且有 \(E_{m2} < E_{m1}\), 在此点有:

    \begin{eqnarray*} a_i + b_i \log(E_{mi}) + c_i \log^2(E_{mi}) & = & a_{i+1} + b_{i+1} \log(E_{mi}) + c_{i+1} \log^2(E_{mi}) \\ b_i + 2 c_i \log(E_{mi}) & = & b_{i+1} + 2 c_{i+1} \log(E_{mi}) \\ \end{eqnarray*}

\(a_i\), \(b_i\), \(c_i\) 为刻度因子,它们与 \(E_{mi}\) 的值由软件刻度得到。如果两条曲线可以拟合到满意的结果,则用不到第三条曲线,此时 \(E_{m2}=0\). SPAN 中会自动计算需要几条拟合曲线。

能量刻度曲线

用于描述指定道能量的函数如下:

\begin{equation*} E_i = a + b \cdot C_i + c \cdot C_i^2 \end{equation*}

此处:

  • \(E_i\) – 第 i 道的能量
  • \(C_i\) – 第 \(C_i\) 道
  • \(a\), \(b\), \(c\) – 刻度系数

如果刻度的点数等于 2, 则因子 c 为 0; 如果则有一个点,则 a, c 均为 0.

实验测量的峰形状

使用已知的标准 γ 源用 γ 射线探测器记录了一系列标准的峰形状。

\begin{equation*} P_i = \frac{A_i - B_i}{S} \end{equation*}

此处:

  • \(P_i\) – 第 i 道的分页因子
  • \(A_i\) – 第 i 道的毛计数
  • \(B_i\) – 第 i 道的本底计数
  • \(S\) – 净峰面积

在 SPAN 中可能使用 64 个点来记录峰形状,这些峰形状存储中库中,能量接近于一条 ROI 的标准峰形状被用来进行非线性最小二乘拟合。标准峰形状还有三个参数:中心道、半高宽 (FWHM) 和高度。标准峰形状和多条 ROI 在比例上不同。需要进行的坐标转换也比较简单。使用实验峰形状可得到更好更快的拟合。

半高宽 (FWHM) 刻度曲线

用于计算半高宽 (FWHM) 的表达式如下:

\begin{equation*} \mbox{FWHM} = \left(a + b\cdot E + c \cdot E^2\right) \cdot \left(1 + d \cdot T_d + e \cdot T_d^2 \right) \end{equation*}

此处:

  • \(E\) – 能量
  • \(T_d\) – 死时间
  • a, b, c, d, e – 刻度因子

高斯形 (Gaussian)

拟合谱的最基本数据表达式。峰在第 i 道的计数 \(F_i\) 由下式给出:

\begin{equation*} F_i = \frac{S}{2.5066 \cdot \sigma} \times e^{-\frac{(i-x_0)^2}{2\sigma^2}} \end{equation*}

此处:

  • \(x_0\) – 拟合峰的中心道
  • \(S\) – 峰面积,与拟合峰的净峰面积相关
  • \(\sigma\) – 峰宽度, \(0.42466\times\mbox{FWHM}\)

SPAN 软件中用户可选择高斯形或实验形来拟合感兴趣的线峰。

线性最小二乘法拟合

许多刻度曲线使用如下形式的多项式:

\begin{equation*} y = a + b \cdot x + c \cdot x^2 \end{equation*}

此处:

  • y – 待计算值 (能量, FWHM 等)
  • x – 需要进行 y 值计算的地方 (道, 能量等)

当得到一组数据时,其因子 a, b, c 可使用最小二乘法拟合得到:

\begin{equation*} \chi^2 = \sum w_i^2 \left( y_i - a - b x_i - c x_i^2 \right)^2 \end{equation*}

此处 \(w_i\) 为权重,通常为 1.0。当 \(\frac{\partial \chi^2}{\partial b}=0\), \(\frac{\partial \chi^2}{\partial c}=0\), \(\frac{\partial\chi^2}{\partial a}=0\), 则可得如下方程组:

\begin{eqnarray*} \sum w_i^2 y_i & = & \sum w_i^2 b + \sum w_i^2 x_i b + \sum w_i^2 x_i^2 b \\ \sum w_i^2 y_i x_i & = & \sum w_i^2 x_i b + \sum w_i^2 x_i^2 b + \sum w_i^2 x_i^3 b \\ \sum w_i^2 y_i x_i^2 & = & \sum w_i^2 x_i^2 b + \sum w_i^2 x_i^3 b + \sum w_i^2 x_i^4 b \\ \end{eqnarray*}

上式可表示为矩阵方程:

\begin{equation*} \sum\limits^2_{j=0} \left(E_{k,j} \cdot f_j\right) = Y_k \end{equation*}

此处:

\begin{eqnarray*} Y_k &= \sum w_i^2 y_i x_i^k &\quad\quad k=0 \sim 2; \\ E_{k,j} &= \sum w_i^2 x_i^{k-j} &\quad\quad j=0 \sim 2; \\ f_j &= a,\, b \quad \mbox{or}\quad c &\quad\quad \\ \end{eqnarray*}

则因子 \(f_j\) 为:

\begin{equation*} f_i = \sum\limits^2_{k=0} \left( E^{-1}_{j,k} \cdot Y_k \right) \end{equation*}

此处 \(E^{-1}_{j,k}\) 为矩阵 \(E_{k,j}\) 的逆。

探测低限 (LLD activity)

对于单峰探测需要达到的最低水平由下式给出:

\begin{equation*} L_d = \frac{4.75 + \sqrt{\sum\limits^{i_0+l}_{i=i_0-l}A_i}}{\gamma \cdot \varepsilon \cdot T_c} \end{equation*}

此处:

  • \(A_i\) – 第 i 道毛计数,由于没有峰,它与本底计数 B_i 相同
  • \(i_0\) – 感兴趣峰的位置
  • \(l = n/2 = 2.5 \cdot \mbox{FWMH}\) – LLD 的积分限
  • \(T_c\) – 采集活时间 (s)
  • γ – 感兴趣 γ 射线的 γ 分支比 (每次退激中出射特定能量的 γ 射线的数量)
  • ε – 感兴趣能点的探测器效率

中子活化分析计算

中子活化分析的绝对方式由下式给出:

\begin{equation*} S = 6.023 \times 10^{23} \cdot c \cdot w \cdot \rho \cdot \left(\Phi\cdot\sigma_{th} + \Phi_e\cdot\sigma_{RI}\right) \cdot \gamma \cdot \varepsilon \cdot T_c \cdot \left(1-e^{-\lambda T_i}\right) \cdot e^{-\lambda T_d} \cdot \left(1-e^{-\lambda T_t}\right) / \left(\lambda \cdot A \cdot T_t \right) \end{equation*}

此处:

  • S – 净峰面积
  • c – 靶元素浓度 (g/g)
  • w – 样品质量 (g)
  • ρ – 靶同位素自然丰度
  • Φ – 热中子能量 (n/s \(\cdot\) cm^2); 热中子平均速度为 2200 m/s (0.0253 eV)
  • σth – 热中子活化截面 ( cm^2 )
  • Φ_e – 镉反应阈能量 (0.5 eV) 以上的中子,通常叫超热或共振中子
  • σRI – 共振积分截面,可由下式定义

    \begin{equation*} \sigma_{RI} = \int\limits^{\infty}_{0.5\mbox{ eV}}\,\sigma(E)\,dE/E \end{equation*}
  • λ – 衰变常数 (s-1)
  • A – 原子质量 (g)
  • T_i – 辐照时间 (s)
  • T_d – 衰变时间 (s)
  • T_c – 采集活时间 (s)
  • γ – 感兴趣 γ 射线的 γ 分支比 (每次退激中出射特定能量的 γ 射线的数量)

参数 S 可由谱分析中得到,参数 Φ , T_i 和 T_d 为样品的辐照与冷却时间,必须在 NAA 计算中给出。采集谱数据时多道分析器会给出参数 T_c 和 T_t . 参数 ε 可由效率刻度得到。其它的参数, λ, ρ, σ, A 和 γ 为保存在数据库中的常量。因此,靶元素浓度 c 可由以前的方程得到。

当然,参数精确时结果亦会正确。但通常获得精确的参数会很困难,比如中了活化截面和探测器效率等,更多的分析使用相对方法。

中子活化分析的相对方法中使用同样的方程式来计算,但计算结果乘以一个修正常数 K 来得到最终结果:

\begin{equation*} C = K \cdot c \end{equation*}

实际修正因子随元素不同亦不同。正确的修正因子可由针对标准参考材料 (SRM) 进行实验得到,并保存于数据库中。当然,修正因子也可由软件中的编辑功能让用户指定。

在 NAA 计算中不同几何位置使用相同的修正因子,因此探测器效率需要针对不同几何位置进行修正,至少相对值需要测量精确。同样,中子能量的相对值也需要进行精确测量。否则,分析结果会不精确。

非线性最小二乘法拟合

在求解多条谱线问题时需要使用非线性最小二乘法拟合,因为峰形状参数在非线性的位置上。形状函数可写为:

\begin{equation*} f(x) = F(h, c, \sigma, x) \end{equation*}

此处:

  • x – 道变量
  • h – 峰高度
  • c – 中心道
  • σ – 峰宽

h, c 和 σ 是待拟合参数,为得到这些参数需要进行如下步骤:

1). 设定参数初始值 h_0, c_0 及 σ_0

2). 计算函数初值及其导数,故函数可写为:

\begin{equation*} f_j(x) = f_{j_0}(x) + \frac{\partial f_{j_0}(x)}{\partial h_j} \Delta h_j + \frac{\partial f_{j_0}(x)}{\partial c_j} \Delta c_j + \frac{\partial f_{j_0}(x)}{\partial \sigma_j} \Delta \sigma_j \end{equation*}

则可得:

\begin{equation*} A_i = B_i + \sum\limits_j f_{j_0}(i) + \sum\limits_j \frac{\partial f_{j_0}(x)}{\partial h_j}\Delta h_j + \sum\limits_j \frac{\partial f_{j_0}(x)}{\partial c_j}\Delta c_j + \sum\limits_j \frac{\partial f_{j_0}(x)}{\partial\sigma_j}\Delta\sigma_j \end{equation*}

此处 \(f_{j_0}(i) = F_j(h_0, c_0, \sigma_0, i)\), 第 j 峰在第 i 道的分布函数初值为 h_0, c_0, σ_0

  • A_i – 第 i 道计数
  • B_i – 第 i 道本底

3). 使用线性最小二乘法拟合得到 Δ h_j, Δ c_j 和 Δσ_j

4). 计算参数初值

\begin{eqnarray*} h^*_0 & = & h_0 + \Delta h \\ c^*_0 & = & c_0 + \Delta c \\ \sigma^*_0 & = & \sigma_0 + \Delta\sigma \\ \end{eqnarray*}

5). 从步骤 2 开始重复,直到 \(\Delta h_j\), \(\Delta c_j\) 和 \(\Delta\sigma_j\) 为 0

多条峰且 σ 变化较小的情况下,则可假设感兴趣区 (ROI) 的多峰有相同的 s 值。如果在 ROI 有 N 条峰,则需要拟合 2N+1 个参数。

寻峰方法

寻峰程序有两条准则:所有的峰都能被找到;光子峰必须能够从本底的波动中区别开来。换句说,不要丢峰,也不要把非峰错误地识别成峰。为达到这两条准则,SPAN 采取两个步骤:

1). 逐道计算来寻找峰,如果 i 道是峰顶端,则需要满足如下条件:

\begin{eqnarray*} A_i > \max(A_{i+1}, A_{i-1}, A_{i-2}, A_{i-3}, A_{i-4}) &\mbox{or}& A_i > \max(A_{i-1}, A_{i+1}, A_{i+2}, A_{i+3}, A_{i+4}) \\ (A_{i+1},A_i,A_{i-1}) > (A_{i+4},A_{i+5},A_{i+6}) + H \cdot k &\mbox{or}& (A_{i+1},A_i,A_{i-1}) > (A_{i-4},A_{i-5},A_{i-6}) + H \cdot k \\ \end{eqnarray*}

此处:

  • max(a,b,c,..) – 即求 a,b,c,… 中的最大值
  • A_i – 第 i 道计数
  • \(k = \sqrt{A_{i-1} + A_i + A_{i+1}}\)
  • H – 寻峰敏感因子。大多数情况下为 3.0, 值越大越敏感,反之则越不敏感。

此步骤可发现所有的单峰和可分开的重峰,但对封闭的重峰中的峰位确定则需要在下一步骤中进行。

2). 重叠的多峰通过检验 FWHM 从而与单峰区别开来。因为重峰比单峰的 FWHM 要大很多。重峰中的弱峰可通过拟合谱并计算余量来识别。

拟合优度 χ

衡量拟合曲线与实验数据符合程序的值为 χ, 对于多峰 ROI 拟合,其值为:

\begin{equation*} \chi = \frac{\sqrt{\sum\limits^r_{k=l}\left(r_k \cdot w_k\right)^2}}{n-m-1} \end{equation*}

此处:

  • n=r-l+1 – ROI 中的总道数
  • r – ROI 的右边界
  • l – ROI 的左边界
  • m – 拟合参数的个数
  • r_k – 第 k 道的余量
  • w_k – 第 k 道的权重

余量

在 ROI 多峰拟合中某道的余量定义如下:

\begin{equation*} r_i = \frac{A_i - B_i - \sum\limits^n_{j=1} F_{ij}}{\sqrt{A_i}} \end{equation*}

此处:

  • A_i – 第 i 道的毛计数
  • B_i – 第 i 道的本底计数
  • Fij – 第 i 道对第 j 峰的拟合计数
  • n – 区域内拟合峰的数目

版权所有 ©2015: Exaos Lee | Date: 2010-12-01 | Generated by Emacs 24.4.1 (Org mode 8.3.2), Validate, 88x31.png

comments powered by Disqus