数据的高维度分析和低维度描述记录-137Cs能谱分析和Eigen

作者: 蒜薹 | 来源:发表于2019-02-16 13:50 被阅读0次

数据的高维度分析和低维度描述记录-137Cs能谱分析和Eigen
数据的高维度分析和低维度描述记录-维度篇
主成份分析
数据建模-维度
第10周降维打击
机器学习day11降维
在线作图|2分钟做Lefse分析
主成份分析算法 PCA
NumPy
维度与能量

数据的高维度分析和低维度描述记录- $^{137}Cs$ 能谱分析和Eigen使用记录

实例- $^{137}Cs$ 能谱

放射源

放射源使用 $^{137}Cs$ ，这是一种同位素，不稳定，容易发生衰变，例如有95%的概率发生 $\beta -$ 衰变，产生 $^{137}Ba^{m}$ 。 $^{137}Ba^{m}$ 是一个所谓的“亚稳态”核，它不稳定，会以85%的概率向自身的基态跃迁。这是一个退激发过程，一个具有动能0.662MeV的光子将会伴随着退激发过程从核中释放出来，我们一般会将这种光子称为 $\gamma$ 射线。

探测器

探测器使用HPGe探测器，它由两个重要部件组成——HPGe晶体和闪烁光转换部件。

HPGE是一种闪烁晶体，由Ge（锗）晶体制作，其中的Ge纯度很高，因此被称作"High-Purity Germanium Detector"（“高纯锗探测器”），简称"HPGe"。锗晶体用来吸收 $\gamma$ 射线的能量，以激发态的形式短暂地存储能量，最终退激发，并将能量以荧光的形式释放。这个转换过程就像晶体在发光一样，由于 $\gamma$ 射线是间断进入锗晶体的，导致晶体是间断着发射荧光的，感觉上像是晶体在“闪烁”。当然，如果 $\gamma$ 射线的强度非常大，在锗晶体完成荧光发射之前，又有新的 $\gamma$ 光子进入了其中，那么锗晶体会继续荧光的发射，整体上就没有了“闪烁”的效果。

闪烁光转换部件，名字叫起来很高大上，其实就是个利用了光电效应的电压表，当年爱因斯坦他们可能就是拿这玩意做的实验。这个东西的主要功能就是把荧光的能量转换为电流的电荷量，中间用了用电子倍增、电磁感应之类的东西，总的来说，就是个光电转换加上一个电压表/电流表。这个东西在测量荧光能量时，会产生误差，而且不小，这也是为啥能谱看上去很圆润的原因。

样本、数据和维度

在本节的实例下，一个样本，表示一个HPGE探测器进行一次测量，测量的时间人为设定好。一个数据，就是一个HPGE进行一次测量的能谱。这里探测器的电子学使用1024道，也就说明一个数据中，包含了1024个数，有1024个维度。

模拟和能谱获取

使用模拟的方法获取能谱。模拟以蒙特卡洛方法为基础，使用建立在这种方法之上的软件Geant4来构建放射源、探测器，并且用它来模拟微观反应的全过程。最后，我来收集 $\gamma$ 射线传递给的能量，并根据这些能量制作能谱。能量的转换过程可以分为两个步骤，第一个步骤是从 $\gamma$ 射线进入锗晶体到锗晶体释放荧光，第二个步骤是荧光进入闪烁光转换部件，转换为电荷量，并以数据形式保存在电脑中。

第一个步骤是 $\gamma$ 光子将自己的动能传递给锗晶体，锗晶体将这些动能以荧光的形式放出，导致锗晶体“闪烁”，这个过程使用Geant4模拟，统计探测器内的总能量沉积即可。

第二个步骤是考虑闪烁光转换部件的测量误差。认为每一次因光子入射而引发的“电压表”记录是一次事件，这次事件的结果是获得一个电荷量，也就是探测器认为的 $\gamma$ 光子传递给它自己的能量。实际上，这个能量是有误差的，一般我们会认为这个误差是以高斯型分布引入的，准确的说，是测量得到的电荷量符合高斯分布，高斯分布的期望是真实的 $\gamma$ 射线传递的能量，标准差是探测器的固有属性，是一个固定的数值。因此，在数学上，这一步的操作就是以上一步得到的能谱为函数，增加一个高斯分布，这种操作属于卷积。

卷积和能谱解释

根据前文的介绍，能谱通过两个步骤制作而成，第一个步骤是统计 $\gamma$ 射线在闪烁体内沉积的能量，以能量沉积谱表征，虽然这里提到了能量沉积谱，但是这个谱并不能在实验中获取，因为它并不是数字信息，不能以数字信息的形式传输到电脑中，更不能被记录在硬盘等数字存储媒介中。为了保存这个能谱，才设计了第二个步骤，简言之就是将第一个步骤中的能量沉积谱转换为数字形式，并且保存在硬盘等数字存储媒介中。第二个步骤是沉积的能量被转换成荧光，荧光接着被转换成电子，电子接着被电场加速，用来影响电极之间的电压，最后电压被电压表测量并以数字信息保存。

虽然步骤二的目标是将能量沉积谱保存在电脑硬盘中，但是步骤二的这一套过程引入了大量的误差，使最终保存的能量不是真实的沉积能量。为了还原真实的能量，需要了解步骤二引入误差的方式。

高斯随机分布，是用来描述步骤二的分布之一，也是最简单的方法。这个方法认为，步骤二得到的能量是一个随机数，这个随机数符合一个高斯分布，这个高斯分布的期望是步骤一的沉积能量，方差与探测器系统有关。因此，步骤一的能量沉积谱中，每一道内的所有计数，也称为“事件”，都符合同一个高斯分布，在步骤二中，它们的能量会按照这个分布被分散到对应的能量区间上。在分散完所有能量道后，就得到了步骤二的结果，这个过程就称为“卷积”。

考察第i个 $\gamma$ 光子入射探测器，也就是第i个事件。该事件中，令光子在闪烁体内沉积的能量为 $EDep_{i}$ ，电子学测量到的能量为 $E_{i}$ 。那么， $E_{i}$ 是一个符合高斯分布的随机数，这个高斯分布的期望 $\mu_{i}$ 为 $EDep_{i}$ ，标准差 $\sigma_{i}$ 与探测器系统有关，并且受到沉积能量的影响。整体上，标准差与探测器系统的关系更大，在探测器系统固定的前提下，不同的能量沉积会造成不同的标准差。在误差允许的前提下，可以认为标准差与能量沉积无关，其实一般都是这么做的。

考察步骤一结束时得到的能谱，这是真正意义上的能量沉积所统计出来的直方图，这些能量是没有被探测器系统“污染”的，没有误差的，不是上文所说的随机数。这个能谱是直方图，就是一系列的区间和计数。设定区间总数量为N，某一个区间使用n表示，注意这个n与上文中的维度编号无关。在本文背景下，我们使用教学常用的1024道电子学系统，那么N为1024。

统计上，认为第 $nBin$ 个区间内的所有事件，都造成了相同的能量沉积，该能量沉积以 $EDep_{nBin}$ 表示，其数值可以以该区间的最小值、最大值、中间值等来代表，本文选择中间值。令 $EDep_{minBin}$ 表示能量沉积直方图中最低道所对应的能量，也就是最小能量沉积所对应的区间； $EDep_{maxBin}$ 表示能量沉积直方图中最高道所对应的能量，也就是最大能量沉积所对应的区间。令 $BinWidth$ 表示一个区间的宽度，根据直方图定义，同一个直方图下，所有区间的宽度相同。
$EDep_{nBin} = EDep_{minBin} + (nBin+0.5)\times BinWidth$

$BinWidth = \frac { EDep_{maxBin} - EDep_{minBin} } {N}$

对于相同区间的事件，也就是能量沉积属于区间编号为 $nBin$ 的事件，它们在步骤一结束时造成的能量沉积在统计上认为相同，我们这里也就认为相同。原因是，由于每一个事件在步骤一的能量沉积数值太多，在模拟中虽然可以全部独立保存，但这样全部独立保存的话，会占用大量的硬盘，而且数据处理起来会消耗大量时间，不实际，因此选择将这些能量沉积数值进行统计，并且以直方图形式保存，也就是上文所谓的“步骤一的能谱”。在这个设定下，一个事件自身的能量沉积数值已经被主动舍弃了，保存下来的只有代表一个区间的能量。也就是说，其实每一个事件的能量沉积数值都可以被保存下来，但是我主动放弃了它们。

考察步骤一的能谱，第 $nBin$ 个区间内的事件，其能量为 $EDep_{nBin}$ 。由于步骤二的测量误差， $EDep_{nBin}$ 将会被施加随机误差，从统计上分析，经过了步骤二的测量，能量值 $E_{nBin}$ 将符合一个高斯分布，这个高斯分布的期望 $\mu_{nBin}$ 为 $EDep_{nBin}$ ，而标准差同样与探测器系统有关，对于同一次测量，认为探测器系统没有改变，由此可认为所有区间的能量 $E_{nBin}$ 所符合的高斯分布的标准差相同。因为第 $nBin$ 个区间内存在不止一个事件，经过了步骤二的误差引入，这些事件的能量将有可能偏离第 $nBin$ 个区间所代表的能量区域。
$E_{nBin} \sim N(EDep_{nBin},\sigma_{nBin})$

$\sigma_{0Bin} = \sigma_{1Bin} = \dots = \sigma_{nBin} = \dots = \sigma_{NBin}$

从统计上分析，步骤一的第 $nBin$ 个区间内的事件，经过步骤二的误差引入后，仍然属于第 $nBin$ 个区间内的概率 $P_{nBin \to nBin}$ ，可以使用高斯分布进行计算。为了更加明确地表示计算，将高斯分布的具体形式引入。
$P_{nBin \to nBin} = \int ^{E_{nBinMax}}_{E_{nBinMin}} {N(EDep_{nBin},\sigma_{nBin})} dE$
为了更加明确地表示计算，将高斯分布的具体形式引入。
$N(EDep_{nBin},\sigma_{nBin}) = \frac {1} {\sqrt{2 \pi \sigma^2}} e^{- \frac {{(E-EDep_{nBin})}^2} {2 \sigma^2_{nBin}}}$
那么， $P_{nBin \to nBin}$ 在真实使用环境下的表达式如下。
$P_{nBin \to nBin} = \int ^{E_{nBinMax}}_{E_{nBinMin}} { \frac {1} {\sqrt{2 \pi \sigma^2}} e^{- \frac {{(E-EDep_{nBin})}^2} {2 \sigma^2_{nBin}}} } dE$

$E_{nBinMax} = EDep_{minBin} + nBin\times BinWidth$

$E_{nBinMax} = EDep_{minBin} + (nBin+1)\times BinWidth$

现在考虑步骤一的第 $nBin$ 个区间内的事件，经过步骤二的误差引入后，属于第 $mBin$ 个区间内的概率 $P_{nBin \to mBin}$ ，其推演与上文相同。
$P_{nBin \to mBin} = \int ^{E_{mBinMax}}_{E_{mBinMin}} { \frac {1} {\sqrt{2 \pi \sigma^2}} e^{- \frac {{(E-EDep_{nBin})}^2} {2 \sigma^2_{mBin}}} } dE$

$E_{mBinMax} = EDep_{minBin} + mBin\times BinWidth$

$E_{mBinMax} = EDep_{minBin} + (mBin+1)\times BinWidth$

上文已经从步骤一的能谱视角出发，介绍了步骤二的误差如何引入。现在从步骤二的能谱视角出发，考虑如何从数学上获得误差干扰后的能谱。

在本文背景下，探测器的电子学获得的能谱是一个离散的数据，对于步骤二能谱的第 $nBin$ 道，其计数，或者说事件，由步骤一的能谱中的所有道贡献而来。由于高斯分布的形状，主要是由第 $nBin$ 道以及附近的若干道贡献。若以高斯分布的 $\pm 3\sigma$ 原则为标准，可以有效的去除无用的区间，提升计算效率。

这里就以 $\pm 3\sigma$ 原则为标准，则步骤二能谱的第 $nBin$ 道由步骤一能谱的第 $nBin$ 道主要构成，同时， $\pm 3\sigma$ 其中的 $(0 \sim 3\sigma)$ 部分，由第 $(n+1)Bin$ 道、第 $(n+2)Bin$ 道至第 $(n+t)Bin$ 道贡献； $\pm 3\sigma$ 其中的 $(-3\sigma \sim 0)$ 部分，由第 $(n-1)Bin$ 道、第 $(n-2)Bin$ 道至第 $(n-t)Bin$ 道贡献。其中，t表示 $(0 \sim 3\sigma)$ 所占有的区间数量，或者是 $(-3\sigma \sim 0)$ 所占有的区间数量，这两种情况所占有的区间数量是一样的，因此统一用t表示。

根据上述分析，步骤二能谱的第 $nBin$ 道计数 $Count_{EnBin}$ ，可以通过上述步骤一的区间来计算。
$Count_{EnBin} = P_{(n-t)Bin \to nBin} \times Count_{EDep(n-t)Bin} + \dots + P_{(n-1)Bin \to nBin} \times Count_{EDep(n-1)Bin} \\ + P_{nBin \to nBin} \times Count_{EDepnBin} \\ + P_{(n+1)Bin \to nBin} \times Count_{EDep(n+1)Bin} + \dots + P_{(n+t)Bin \to nBin} \times Count_{EDep(n+t)Bin}$

$Count_{EnBin} = \sum _{i=n-t}^{n+t} {P_{nBin \to iBin} \times Count_{EDepiBin}}$

以矩阵形式表达更加简明。
$Count_{EnBin} = \begin{bmatrix} 0 & \dots & 0 & P_{(n-t)Bin \to nBin} & \dots & P_{nBin \to nBin} & \dots & P_{(n+t)Bin \to nBin} & 0 & \dots & 0 \end{bmatrix} \begin{bmatrix} 0 \\ \vdots \\ 0 \\Count_{EDep(n-t)Bin} \\ \vdots \\ Count_{EDepnBin} \\ \vdots \\Count_{EDep(n+t)Bin} \\ 0 \\ \vdots \\ 0 \end{bmatrix}$

img

能谱数据处理

中心化

中心化是针对样本矩阵而言的，针对每一个维度，将数据剪去期望。

在本文背景下，样本数据的每一列代表一个维度，每一行代表一个样本。在求取期望时，需要知道行数和列数，在Eigen3下，获取行数和列数的方法如下。

int numRows = SampleMatrix_.rows();
int numCols = SampleMatrix_.cols();

协方差矩阵

协方差矩阵由中心化的样本矩阵计算得到。

协方差矩阵的图形显示可以通过“热力图”来完成。“热力图”实际上就是一个二维直方图，使用色彩来表征每一个区间(bin)的计数值。Python下的Seaborn库中有热力图的支持，使用起来非常简单，可参考[Jianshu Matplotlib, 2017]，[Seanborn, Heatmap, 2019]。Seaborn是以Matplotlib为基础的可视化库，比Matplotlib使用更加简洁方便。

特征值和特征向量

Eigen库中的类SelfAdjointEigenSolver专门负责特征值和特征向量求解，具体参考[Eigen-EigenDeconposition, 2019]。

并不是任何矩阵都可以进行特征分解的，在本文应用背景下，由样本数据组成的矩阵只有很小概率可以进行特征分解，而且，无论是在PCA还是因子分析中，样本矩阵都不是用来进行特征分解的对象。

在PCA中，进行特征分解的矩阵是样本矩阵的协方差矩阵。协方差矩阵是一个对称阵，而且所有元素都是实数，一般称为实对称矩阵，这种矩阵的特点之一就是可以进行特征分解。

MatrixXd ones = MatrixXd::Ones(3,3);
SelfAdjointEigenSolver<MatrixXd> es(ones);
cout << "The first eigenvector of the 3x3 matrix of ones is:" 
     << endl << es.eigenvectors().col(1) << endl;

参考文献

网站

[Eigen-EigenDeconposition, 2019] Eigen官方说明文档，特征分解类介绍 http://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html#title17

[Jianshu Matplotlib, 2017] 用python-pandas作图矩阵 https://www.jianshu.com/p/47b54eb35eed

[Seanborn, Heatmap, 2019] Seaborn API for heatmap, http://seaborn.pydata.org/generated/seaborn.heatmap.html#seaborn.heatmap

论文

[Miller, E. A., 2015]: Miller, E. A. , Robinson, S. M. , Anderson, K. K. , Mccall, J. D. , Prinke, A. M. , & Webster, J. B. , et al. (2015). Adaptively reevaluated bayesian localization (arbl): a novel technique for radiological source localization. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 784, 332-338.

[Penny, R. D., 2015]: Penny, R. D. , Crowley, T. M. , Gardner, B. M. , Mandell, M. J. , Guo, Y. , & Haas, E. B. , et al. (2015). Improved radiological/nuclear source localization in variable norm background: an mlem approach with segmentation data. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 784, 319-325.

[Jianshu Matplotlib, 2017]:

数据的高维度分析和低维度描述记录-137Cs能谱分析和Eigen
数据的高维度分析和低维度描述记录-能谱分析和Eigen使用记录实例-能谱放射源放射源使用，这是一种同位素，不...
数据的高维度分析和低维度描述记录-维度篇
数据的高维度分析和低维度描述记录-维度篇背景和目标介绍本文是我近期展开的一个研究——放射源定位——方向的思考和...
主成份分析
概念主成份分析（PCA）是一种维度归约的技术。（维度规约：把高维数据转换成低维数据的。）涉及的其他概念数据标...
数据建模-维度
维度的基本概念维度是维度建模的基础和灵魂。在维度建模中，将度量称为“事实”将环境描述为“维度”，维度是用于分析事...
第10周降维打击
降维打击原意指的是不同维度的竞争者，高维度生物毁灭低维度生物，却与高维度生物无关。 1. 竞争者假设A和B在竞争...
机器学习day11降维
降维用一个低维度的向量表示原来高维度的特征，避免维度灾难。降维方法主成分分析线性判别分析等距映射局部线...
在线作图|2分钟做Lefse分析
Lefse分析 Lefse（LDA Effect Size）分析是一种用于发现和解释高维度数据的基因、通路和分类单...
主成份分析算法 PCA
PCA 算法主要是把高维度的数据降为低维度数据。典型地应用包括数据压缩和数据可视化。本文介绍 PCA 算法及其典型...
NumPy
一、数据的维度数据维度的python表示一维数据：列表和集合类型二维数据：列表类型多维数据：列表类型高维数...
维度与能量
维度与能量的关系结论越高维度需要消耗的能量越高低维度具有低能耗的特性成本低廉等特性低维度进入高纬度？高纬...

数据的高维度分析和低维度描述记录-137Cs能谱分析和Eigen

数据的高维度分析和低维度描述记录- $^{137}Cs$ 能谱分析和Eigen使用记录

实例- $^{137}Cs$ 能谱

放射源