基于压缩感知理论的层析γ扫描技术的透射图像重建方法与流程

基于压缩感知理论的层析
γ
扫描技术的透射图像重建方法
技术领域
1.本发明属于核废物γ射线无损检测技术，具体涉及一种基于压缩感知理论的层析γ扫描技术的透射图像重建方法。


背景技术：

2.对核废物进行基于γ射线的无损检测(non-destructive assay,nda)技术主要包括分段γ扫描(segmented gamma scanning,sgs)技术和层析γ扫描(tomographic gamma scanning,tgs)技术。tgs技术能够精确检测非均匀分布的核废物，可精确定位核废物中的γ放射性核素并识别其种类、含量。tgs技术将装载核废物的容器划分为多层，每层进一步划分为若干个体素(立方体元)并假定每个体素内介质和放射性均匀分布。tgs技术对核废物的扫描测量分为透射测量和发射测量，通过透射测量重建得到各层的二维透射重建图像(核废物的线性衰减系数(μ,cm-1
)分布图)，整合整桶所有层的透射重建结果后，为发射测量提供体素介质对体素源的衰减校正以获得最终的三维发射重建图像(核废物的放射性活度(s,bq)分布图)。作为一种实用型工业技术，较长的测量时间阻碍了tgs技术的大规模商业应用，可通过连续扫描模式降低其测量时间。在连续扫描模式下，信噪比较低或者检测的核废物密度过高等导致所测的γ射线数据统计性不足，采样数据不能满足传统的nyquist-shannon(whittaker nyquist kotelnikov shannon)采样定理，需要对tgs技术进行欠采样数据下的不适定(ill-posed)数学反问题的求解。2006年压缩感知(compressed sensing,cs)理论被candes、donoho等人提出并于同年被sidky等人成功应用于计算机断层扫描(computed tomography,ct)技术的图像重建中，后者证实该理论可基于图像的稀疏性通过极少的测量数据获得质量较优的重建图像，为这一问题的解决提供了新的思路。sidky等人在代数重建算法(algebraic reconstruction technique,art)迭代过程中引入图像总变分/总变差(total variation,tv)最小化(minimization)约束，该算法即为经典的art-tv算法，又被称为tv-art算法或art-tvm算法。
3.cs理论成功应用于图像重建的前提条件是图像具有稀疏性。ct图像中，人体内的相同组织具有近似的衰减系数，其ct图像灰度值近似分片光滑，利用这一特点将ct图像进行离散梯度变换(discrete gradient transform,dgt)得到的图像梯度通常是一个稀疏图像，对图像梯度的稀疏性约束就形成了常用的图像tv最小化约束。将cs理论应用于tgs技术的二维透射图像重建时，也需要tgs透射图像具有稀疏性。图像的稀疏性与图像的稀疏度有关。图像的稀疏度是指图像中的非零元素的个数与元素总数的比值。图像中的非零元素越少，图像的稀疏度数值越小，图像的稀疏性越高，图像就能以较大概率精确重建且所需的测量数据的个数越少。由于核废物内介质分布未知，为将cs理论成功应用于tgs二维透射图像的重建，需将tgs图像进行变换以稀疏化处理。因此，基于tgs技术的每个体素内的介质均匀分布的假定，对tgs二维透射图像进行提高图像分辨率的设置后，每个体素被进一步划分为多个更小的体素，则提高图像分辨率后的每个更小的体素内均按照假定与原体素具有相同的衰减系数，然后对其图像梯度进行稀疏性变换，满足了cs理论的图像稀疏性的应用要求。
基于cs理论，将图形tv最小化约束应用于图像重建中时，图像tv最小化需要与迭代算法相互结合才能进行应用。迭代算法是把扫描测量的数据以离散化的数学模式用线性方程组表示，然后迭代计算求解相应图像。对于tgs技术，比较有代表性的迭代算法有代数重建算法(algebraic reconstruction technique,art)和极大似然期望最大化(maximum likelihood expectation maximization,mlem)算法。联合代数重建算法(simultaneous algebraic reconstruction technique,sart)和联合迭代重建算法(simultaneous iterative reconstruction technique,sirt)都是art算法的改进算法。有序子集期望最大化(ordered subsets expectation maximization,osem)算法则是mlem算法的改进算法。art算法不属于但mlem算法属于统计迭代算法，统计迭代算法考虑了投影数据的统计特性，即投影数据满足泊松分布，但在实际应用中，投影采集到的测量数据可能并不严格遵守高斯或泊松分布，尤其是在欠采样数据的情况下，因此可选用art算法或其改进算法作为与图像tv最小化相结合的迭代算法，得到art-tv算法，sart-tv算法或sirt-tv算法等。


技术实现要素：

4.本发明的目的是针对现有技术中存在的问题，提供一种基于压缩感知理论的层析γ扫描技术的提高图像分辨率的透射图像重建方法，以提高tgs技术的透射重建图像的质量。
5.本发明的技术方案如下：一种基于压缩感知理论的层析γ扫描技术的透射图像重建方法，包括如下步骤：
6.(1)计算透射源发射的特征γ射线穿过核废物桶中一层体素时的透射径迹长度对应的矩阵t，根据透射重建原理获得透射测量数据并处理后得到透射源发射的特征γ射线穿过该层体素的透射率的负自然对数对应的矩阵p；
7.(2)预设art算法的迭代次数iter
art
和迭代松弛因子λ
art
，对预设的透射重建图像赋初值μ0，采用art算法按照预设的迭代次数进行迭代计算，得到透射重建图像μj；
8.(3)预设图像tv迭代的迭代次数iter
tv
和迭代松弛因子λ
tv
，针对核废物桶二维透射图像全部体素，对每一次art重建后得到的μj进行图像tv调整，使迭代重建图像的tv最小化；
9.(4)将图像tv最小化后处理的图像提供给art算法作为迭代的初值，重复上述步骤(2)和(3)，直到满足设定的收敛条件(如：相隔10000次总迭代的数据差值绝对值的最大值＜1e-7)后自动停止迭代。
10.进一步，如上所述的基于压缩感知理论的层析γ扫描技术的透射图像重建方法，步骤(2)中art算法的迭代公式为：
[0011][0012]
式中，
[0013]
j为迭代次数，j≤iter
art
，其中iter
art
为预设的art算法的迭代次数，
[0014]
和分别为第m行，第n列体素的线性衰减系数的最新迭代值和当前迭代
值，
[0015]
ti为行向量，是透射径迹长度矩阵t＝[t1,t2,
…
,ti]
t
的第i行，代表着第i条射线，
[0016]
pi为透射率的负自然对数矩阵p的第i行数据，
[0017]
为考虑第i条射线影响后的当前估计值。
[0018]
进一步，如上所述的基于压缩感知理论的层析γ扫描技术的透射图像重建方法，步骤(2)中art算法的迭代计算约束条件为各体素的μ
m,n
的非负性。
[0019]
进一步，如上所述的基于压缩感知理论的层析γ扫描技术的透射图像重建方法，步骤(3)中，在各层桶外添加体素以保证核废物桶二维透射图像全部体素参与图像tv最小化迭代，所添加的体素包括各层桶外四个边的体素以及四个角的体素。
[0020]
更进一步，对于所添加的各层桶外四个边的体素，引用neumann边界条件进行赋值，其值为与其同层相邻的桶内的体素的数值；
[0021]
对于所添加的各层桶外四个角的体素，其值为与其同层最相近的靠近桶边缘的两个体素的数值的平均值。
[0022]
进一步，如上所述的基于压缩感知理论的层析γ扫描技术的透射图像重建方法，步骤(3)中，图像tv最小化的公式如下：
[0023][0024]
式中，
[0025]
d为图像tv下降的梯度的步长，为透射图像tv的导数。
[0026]
更进一步，d为固定值，或者也可以对其进行改动以使其自适应下降；可采用图像tv梯度自适应下降法，设置d为经过图像tv约束后得到的图像与未经图像tv约束得到的图像的差值，使其自适应下降。
[0027]
本发明的有益效果如下：本发明通过对tgs技术的二维透射图像进行提高分辨率的设置，使提高分辨率后的tgs二维透射图像梯度具有稀疏性，进而基于cs理论，采用art-tv算法作为提高分辨率后的tgs二维透射图像的重建算法，并引入neumann边界条件和其他设置使tgs二维透射图像中的所有体素参与图像tv最小化迭代过程，从而提升tgs技术的透射重建图像的质量。
附图说明
[0028]
图1为核废物桶单层88(对应于10
×
10划分方式)个体素的分布示意图；
[0029]
图2为二维图像tv最小化时单个体素μ
m,n
所涉及的周围体素分布示意图；
[0030]
图3为整层344(对应于20
×
20划分方式)个体素在透射图像tv最小化时的分布示意图；
[0031]
图4为整层344(对应于20
×
20划分方式)个体素在透射图像tv最小化时全部进行迭代时的体素分布示意图；
[0032]
图5为整层所有体素参与透射图像tv最小化时桶外添加体素说明截取图(取自图4右上角)；
[0033]
图6为整层756(对应于30
×
30划分方式)个体素在透射图像tv最小化时全部进行迭代时的体素分布示意图；
[0034]
图7为整层1324(对应于40
×
40划分方式)个体素在透射图像tv最小化时全部进行迭代时的体素分布示意图。
具体实施方式
[0035]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0036]
对于tgs二维透射重建图像问题，运用cs理论时对应的图像重建问题就是寻求下述最稀疏解的问题：
[0037]
min||μ||
tv
,s.t.t
·
μ＝p,μ
m,n
≥0
[0038]
式中，t为透射源发射的特征γ射线穿过某一层体素时的透射径迹长度对应的矩阵，其实质上是一个大型的稀疏矩阵；μ为待求的该层体素对应于该特征γ射线的线性衰减系数矩阵；p为透射源发射的该特征γ射线穿透该层体素的透射率的负自然对数对应的矩阵。
[0039]
图像梯度算子一般用符合表示，对于tgs二维透射图像的各个体素(m,n)，其图像梯度与其周围的体素有关，其各维的图像梯度为：
[0040][0041]
即，tgs二维透射图像中体素对应的图像梯度值为：
[0042][0043]
式中，μ
m,n
为二维透射重建图像中第m行，第n列体素对应于某一特征能量的γ射线的线性衰减系数(cm-1
)。
[0044]
其二维透射图像tv为：
[0045][0046]
实际上，为方便进行求导计算，通常采用l2范数近似代替上述的l1范数，此外为避免图像tv求导后分母为零或无穷大，需引入一个较小的正数ε，一般ε＝10-8
。即
[0047][0048]
上式对某一体素μ
m,n
进行求导，其对应偏导数为：
[0049][0050]
以采用tgs技术检测208l的标准钢桶为例，若采用目前国际上流行的基于直角坐标系下的体素划分方式，每层按照桶内径划分为10
×
10个体素，准直器开口与体素大小相当，不考虑每层四个边角处各位于桶外的3个体素，即每层共88个体素，包括60个规则形状的立方体体素和位于边界处的28个不规则形状体素，每层体素分布如图1所示。
[0051]
此处为演示cs理论在提升图像分辨率后的tgs二维透射重建图像的应用情况，以每层344(对应于20
×
20划分方式)个体素为例进行说明。由上述公式可知，在进行二维透射图像tv最小化时，单个体素μ
m,n
的处理需涉及如图2所示的其周围的6个体素。因此，整层344个体素在tgs二维透射图像tv最小化过程中，共涉及276个体素，如图3所示。
[0052]
需要注意的是，上述tgs二维透射图像tv最小化过程中只涉及非边界体素，且这些非边界体素不与边界接触。该种图像tv最小化过程不涉及边界体素的情况实际上满足ct技术中需要保持图像边缘的要求，但是对于tgs技术而言并没有保持图像边缘的要求，而只需要在透射图像重建中考虑桶内核废物的衰减特性的精确确定。因此，在tgs二维透射图像tv最小化过程中可通过考虑整层的所有体素以提高所获得的透射图像质量，此种情况下的图像tv最小化时对应的体素分布如图4所示。
[0053]
以图4右上角为例进行说明，见图5，有：
[0054]
(1)对于图像中桶外标注的体素，如图中体素μ
0,15
和μ
6,21
，这些体素引用neumann边界条件(neumann boundary conditions)进行赋值，有：
[0055][0056]
即有：
[0057][0058]
(2)对于图像中桶外标注
①
的体素，neumann边界条件不再适用，定义其值为与其最相近的最靠近桶边缘的两个体素的数值的平均，其值需在标注的体素和桶内体素的基础上得到，即有：
[0059][0060]
(3)对于图像中桶外标注
②
的体素，neumann边界条件不再适用，定义其值为与其最相近的最靠近桶边缘的两个体素的数值的平均，其值需在桶外标注
①
，或者
①
和的体素的基础上得到，即有：
[0061][0062]
(4)对于图像中桶外标注
③
的体素，neumann边界条件不再适用，定义其值为与其最相近的最靠近桶边缘的两个体素的数值的平均，其值需在桶外标注
②
，或者
②
和
①
的体素的基础上得到，即有：
[0063][0064]
上述以基于直角坐标系下的单层的20
×
20划分方式为例。实际上，为进一步提高tgs技术的二维透射重建图像的质量，在应用cs理论时也可进一步提升透射图像的分辨率，如可进一步将单层划分为30
×
30，40
×
40等多个体素以进行基于cs理论下的tgs二维透射图像重建，其对应的整层的所有体素在二维透射图像tv最小化中全部进行迭代时的体素分布分别如图6和图7所示。
[0065]
基于压缩感知理论，采用art-tv算法进行tgs二维透射图像重建时，其一次总迭代的循环步骤可分为两层：(1)外层迭代循环，采用art算法进行迭代重建；(2)内层迭代循环，实现tgs二维透射图像tv最小化。然后将图像tv最小化后处理的图像赋值给art算法作为迭代初值，进行下一次迭代。交替进行art算法迭代和tgs二维透射重建图像tv最小化，直到满足设定的收敛条件后自动停止迭代。其算法的具体实施过程如下：
[0066]
min||μ||
tv
,s.t.t
·
μ＝p,μ
m,n
≥0
[0067]
(1)计算透射径迹长度矩阵t，根据透射重建原理获得透射测量数据并处理后得到
p矩阵。t为透射源发射的特征γ射线穿过某一层体素时的透射径迹长度对应的矩阵，其实质上是一个大型的稀疏矩阵；μ为待求的该层体素对应于该特征γ射线的线性衰减系数矩阵；p为透射源发射的该特征γ射线穿透该层体素的透射率的负自然对数对应的矩阵。
[0068]
(2)预设art算法的迭代次数iter
art
和迭代松弛因子λ
art
，对预设的透射重建图像赋初值μ0，采用art算法按照预设art算法的迭代次数进行迭代计算，得到透射重建图像μj，约束条件为各体素的μ
m,n
的非负性。对应的art算法的迭代公式为：
[0069][0070]
式中，j为迭代次数(j≤iter
art
，其中iter
art
为预设的art算法的迭代次数)；和分别为第m行，第n列体素对应于该特征γ射线的线性衰减系数的最新迭代值和当前迭代值；ti是行向量，是透射径迹长度矩阵t＝[t1,t2,
…
,ti]
t
的第i行，代表着第i条射线；pi是上述的透射率的负自然对数矩阵p的第i行数据，为考虑第i条射线影响后的当前估计值。
[0071]
(3)预设图像tv迭代的迭代次数iter
tv
和迭代松弛因子λ
tv
，对每一次art重建后得到的μj利用梯度下降法进行tv调整，使迭代重建图像的tv最小化，对应公式如下：
[0072][0073]
式中，d为图像tv下降的梯度的步长，可使其固定不变或也可对其进行改动以使其自适应下降，为透射图像tv的导数。可采用图像tv梯度自适应下降法，设置d为经过图像tv约束后得到的图像与未经图像tv约束得到的图像的差值，使其自适应下降。
[0074]
(4)将图像tv最小化后处理的图像提供给art算法作为迭代的初值，重复上述步骤(2)和(3)，直到满足设定的收敛条件(相隔10000次总迭代的数据差值绝对值的最大值＜1e-7)后自动停止迭代。
[0075]
本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本发明的其它实施方案。本技术旨在涵盖本发明的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本发明的一般性原理并包括本发明未公开的本技术领域中的公知常识或惯用技术手段。这样，倘若对本发明的这些变型、用途适应性变化属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改变型和用途适应性变化在内。
[0076]
上述实施方式只是对本发明的举例说明，本发明也可以以其它的特定方式或其它的特定形式实施，而不偏离本发明的要旨或本质特征。因此，描述的实施方式从任何方面来看均应视为说明性而非限定性的。本发明的范围应由附加的权利要求说明，任何与权利要求的意图和范围等效的变化也应包含在本发明的范围内。