CN103207946A - 基于截断奇异值和全变分的闪光照相客体正则化重建方法 - Google Patents

Abstract

本发明公开了一种基于截断奇异值和全变分的闪光照相客体正则化重建方法，包括：（1）通过正向光子输运程序，获取闪光照相成像底片信息；（2）通过截断奇异值方法获得闪光照相客体线性系数空间分布和几何参数；（3）通过全变分正则化方法获得闪光照相客体线性系数空间分布和几何参数。该方法能有效重构闪光照相客体的线性吸收系数，具有较高的计算精度，同时在客体内部材料交界处具有清晰的分辨，具有较高的边缘保真能力，满足工程要求。

Description

基于截断奇异值和全变分的闪光照相客体正则化重建方法

技术领域

本发明属于光子成像技术领域，涉及一种闪光照相客体重建方法，尤其是一种基于截断奇异值和全变分的闪光照相客体正则化重建方法。

背景技术

高能闪光照相利用高能射线穿透客体成像。在高能闪光照相系统中，电子在直线感应加速器加速后，撞击靶件，产生高能（MeV）光子，利用高能光子照射客体，高能光子在客体中进行输运过程，并与之发生相互作用，如康普顿散射，电子对效应和光电效应等，最后穿透客体，并在底片成像。散射照射量通过一定的设备扣除，通过研究透射粒子在底片的成像信息来反推客体的物理性质和几何性质。故能否从底片透射照射量信息建立合适的重建模型对精确反演客体的物理参数和几何参数影响很大。基于此，对高能闪光照相系统重建模型进行研究，并对模型的稳定性进行研究，对于提高重构客体质量、获得精确的线性吸收系数的空间分布和客体的几何结构分布具有重要意义。

目前国内外开展了一些闪光照相客体重建研究，其中包括解析重建和迭代重建。解析重建包括的方法有滤波放投影等技术，它先对投影数据进行滤波，再把滤波后的投影数据进行反投影计算。这类算法具有分辨率高的特点，但是它对数据完备性具有严格要求，这就意味着需要进行较长时间和全角度范围的检测，这带来了昂贵的检测成本。同时，滤波反投影算法对数据的噪声敏感性较大，如果探测数据存在散射噪声，重建结果严重偏离真值。迭代重建算法有明确的几何意义和物理意义，而且重建算法简单，但是迭代重建算法计算量大，计算速度相对较慢，并且其对初值的选择也比较敏感，较差的初值选取严重影响重建结果。在迭代算法中，最为熟悉的是代数重构算法。

基于正则化的重建是今年来研究的热点，这类算法主要有Tikhonov正则化算法，Landweber迭代正则化算法。Tikhonov正则化提出较早，在图像处理和信号处理领域运用较多。在闪光照相系统中，采用锥形束投影，一次采集可以获取数目较大的采集量。传统的最小二乘算法可以保证重建系统获得物理解，然而在闪光照相系统中，光子与照相器件发生散射，这个底片带来严重散射噪声，不仅减低了图像的对比度，而且在反演过程中重建结果的严重振荡。Tikhonov在最小二乘的基础上加上了带有正则化参数的惩罚正则化项，虽然可以抑制噪声，但是惩罚项是基于二范数的，它严重平滑了重建客体的边界信息，给边缘带来模糊。Landweber迭代正则化不像Tikhonov一样采用一些复杂的准则确定正则化参数，而是采用迭代步数的倒数作为正则化参数，这大大减少了系统的计算量，然而Landweber迭代正则化算法出现半收敛，在迭代次数达到一定数目时，系统计算收敛，随着迭代次数的增加，系统又趋于发散。Landweber迭代正则化算法中的参数选取过于依赖经验。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种基于截断奇异值和全变分的闪光照相客体正则化重建方法，其一方面能够获得精确的重构结果，消除边缘模糊，另一方面能够抑制闪光照相系统散射噪声对重建结果的放大作用，同时正则化参数通过先验选取，不依赖经验。

本发明的目的是通过以下技术方案来解决的：

该种基于截断奇异值和全变分的闪光照相客体正则化重建方法，包括以下步骤：

（1）通过正向光子输运程序，获取闪光照相成像底片信息；

（2）通过截断奇异值方法获得闪光照相客体线性系数空间分布和几何参数；

（3）通过全变分正则化方法获得闪光照相客体线性系数空间分布和几何参数。

进一步的，步骤1）具体包括以下步骤：

1）从成像底片上选取获得各个像素点的中心坐标；

2）光子输运模型计算得到成像底片像素点透射照射量信息；

3）从闪光成像底片透射照射量，利用闪光照相截断奇异值方法得到闪光照相客体吸收系数空间分布；

4）从闪光成像底片透射照射量，利用闪光照相全变差正则化方法得到闪光照相客体吸收系数空间分布。

上述的闪光照相截断奇异值方法为：从投影矩阵直接出发，获得投影矩阵的奇异值分解，并将投影矩阵数值小的奇异值截断，保留数值大的奇异值；实现对含散射噪声重构时误差放大的抑制。

上述光子输运模型为光子输运正向模型。

上述的闪光照相全变差正则化方法为：从投影的最小二乘方程出发，添加L1范数的惩罚项，并予以正则参数平衡最小二乘项和正则项的平衡关系；在获得高精度的重构结果的同时，具有边缘保真能力。闪光照相正则参数的确立方法为：利用数学准则，最优化确定正则参数；所述数学准则包括广义交叉验证原理，Ｌ曲线原理，准最优化。

与现有技术相比，本发明有如下突出优点：

1.针对现有闪光照相技术，提出了截断奇异值正则化和全变差正则化的重构技术。克服了闪光照相在重构过程中出现的不适定现象，抑制了散射噪声对重建结果误差的放大作用。

2.正则参数通过数学准则先验选取，减少对经验关系的依赖，节约了累计实验所需的大量成本，拓展了重建系统的应用范围。

3.相对代数迭代算法，计算量减低，有效的提高了计算效率。

4.相对传统的Tikhonov方法，本重建系统对客体几何边缘具有较高的保真能力。为闪光照相客体性质诊断提供依据。

附图说明

图1为闪光照相系统示意图；

图2为含散射修正图像重建系统流程图。

具体实施方式

本发明用于提供一种闪光照相中，对客体线性吸收系数空间分布和几何参数的空间分布进行重构的截断奇异值方法和全变分正则化方法。重建系统包括正向光子输运计算程序，获取成像底片透射照射量信息。从底片透射照射量出发，进行反演重构闪光照相客体的物理参数和几何参数。底片记录了大量的客体信息，重构问题是一个不适定的问题，需要诉诸正则化方式进行求解。截断奇异值方法将闪光照相投影矩阵进行奇异值分解之后数值较小的奇异值截断，从而有效的抑制系统对底片散射噪声的波动放大作用。针对系统线性吸收系数空间分布不连续的特点，在最小二乘项的基础上增加了全变分正则化项，并且加以适当的正则参数以平衡残差项和惩罚项的比例。

本发明是基于图1的闪光照相系统进行的，该基于截断奇异值和全变分的闪光照相客体正则化重建方法，具体包括以下步骤：

（1）通过正向光子输运程序，获取闪光照相成像底片信息；具体包括以下步骤：

1）从成像底片上选取获得各个像素点的中心坐标；

2）光子输运正向模型计算得到成像底片像素点透射照射量信息；

3）从闪光成像底片透射照射量，利用闪光照相截断奇异值方法得到闪光照相客体吸收系数空间分布；所述的闪光照相截断奇异值方法为：从投影矩阵直接出发，获得投影矩阵的奇异值分解，并将投影矩阵数值小的奇异值截断，保留数值大的奇异值；实现对含散射噪声重构时误差放大的抑制。

4）从闪光成像底片透射照射量，利用闪光照相全变差正则化方法得到闪光照相客体吸收系数空间分布。所述的闪光照相全变差正则化方法为：从投影的最小二乘方程出发，添加L1范数的惩罚项，并予以正则参数平衡最小二乘项和正则项的平衡关系；在获得高精度的重构结果的同时，具有边缘保真能力。正则参数的确立方法为：利用数学准则，最优化确定正则参数；所述数学准则包括广义交叉验证原理，Ｌ曲线原理，准最优化。

本发明的最佳实施例中，采用的光子输运模型为光子输运正向模型。

（2）通过截断奇异值方法获得闪光照相客体线性系数空间分布和几何参数；

（3）通过全变分正则化方法获得闪光照相客体线性系数空间分布和几何参数。

要实现以上方法，本发明包括以下模块：

闪光照相透射照射计算模块，基于快速光子输运过程获取闪光照相底片透射照射量强度分布图；

重建模块，对上述透射照射量分布进行闪光照相图像重建，针对重建问题的不适定性建立截断奇异值正则化模型和全变差模型，重构出闪光照相客体材料的线性吸收系数信息或者密度信息。

其中透射照射量模块采用快速光子输运计算解析方法，在光子输运过程中，假设光子在于客体发生碰撞，散射光子主要来自初级散射光子，且对于几何结构确定散射光子份额基本不变。光子一经散射，即被认为消失。所以在求解透射照射量的过程中采取精确的光子输运程序。在获取精确的透射照射量之后，完成了对成像底片的模拟。通过底片信息，建立重构模型，采用截断奇异值正则化和全变差正则化方法,采取计算准则得到出的先验正则化参数计算客体的线性吸收系数。

本发明中，截断奇异值的思想是将闪光照相投影方程中投影矩阵进行奇异值分解，将数值较小的奇异值截断，保留数值较大的奇异值，抑制投影矩阵较小奇异值对重构线吸收系数带来的波动放大。全变分正则化方法的思想是在原投影方程的最小二乘项的基础上，加上全变分正则化项。利用全变分正则化项具有较好的处理不连续线吸收系数的性质，对闪光照相客体内部相邻两种材料的分界做出明显界定。

图2表示针对闪光照相系统，不含散射情况下，重建系统截断奇异值和全变分正则化重建的方法。包括以下的步骤：

1）针对已知客体信息，由快速光子输运计算，计算得到模拟成像底片信息中的透射量。

2）由上述得到的透射照射量，进行闪光成像图像重建模块，完成客体的线性吸收系数和几何边界正则化重建。

1.针对上述步骤中的光子输运计算模块，其具体步骤如下：

高能光源发射出粒子，粒子在介质中衰减，衰减关系满足朗伯定律，

I＝I₀exp(-∫Σldl) （1）（1）式可以整理为，

$I=I_0\exp (-\underset{k}{\mathrm{\Σ}}{\mathrm{\Σ}}_k{d}_{i,k})---\left(2\right)$

从光源到成像底片像素点连线上的每条线都满足（2）式，把（2）式写成矩阵形式：

$I_0\exp \{-\left[\begin{array}{ccccc}{d}_{\mathrm{1,1}}& d_{\mathrm{1,2}}& d_{\mathrm{1,3}}& \·\·\·& d_{1,\mathrm{ne}}\\ d_{\mathrm{2,1}}& d_{\mathrm{2,2}}& d_{\mathrm{2,3}}& \·\·\·& d_{2,\mathrm{ne}}\\ \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·\\ d_{\mathrm{number},1}& d_{\mathrm{number},2}& d_{\mathrm{number},3}& \·\·\·& d_{\mathrm{number},\mathrm{ne}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\Σ}}_1\\ {\mathrm{\Σ}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\Σ}}_{\mathrm{ne}}\end{array}\right]\}=\left[\begin{array}{c}{I}_1\\ I_2\\ \·\\ \·\\ \·\\ I_{\mathrm{ne}}\end{array}\right]---\left(3\right)$

其中I₀、I分别为光子束穿过客体前后的光子强度；Σ_k、d_i,k分别为第k层客体的线性吸收系数和射线穿过该客体层的几何距离；number为底片像素点数;ne为客体材料分区数目。

2.在获得底片透射照射量的信息前提下，实现重建。在闪光照相中，采取锥形束投影，底片信息足够多的情况下，

$\underset{k}{\mathrm{\Σ}}{d}_{i,k}{\mathrm{\Σ}}_k=\ln \left(\frac{I_0}{I}\right)---\left(4\right)$

方程（4）可以简写为算子方程

L^m×nΣ^n×1＝B^m×1  （5）

其中m＞n，这是一个超定方程组，由于其不适定性，通常的线性代数手段是无法进行求解的，需要借助于最小二乘法，

$\min {\left|\right|\mathrm{L\Σ}-B\left|\right|}_2^2---\left(6\right)$

若L的奇异值分解为

$L=U\left[\begin{array}{cc}\mathrm{\Φ}& 0\\ 0& 0\end{array}\right]V^T---\left(7\right)$

系数矩阵L的广义逆为

L⁺＝V₁Φ^-1U₁ ^T  （8）

方程族的最小二乘解为，

${\mathrm{\&Sigma;}}_{\mathrm{ls}}=L^{+}B=V_1{\mathrm{\&Phi;}}^{-1}{U_1}^{T}B=\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}\frac{<u_i,b>}{{\mathrm{\&phi;}}_i}{v}_i---\left(9\right)$

其中φ_i为系数矩阵的奇异值，随着下标的增加，系统的奇异值迅速减小。从方程式（9）可以看出，右端项B存在噪声，即底片中存在成像系统的散射照射量，噪声将被放大，使得重建结果偏离真实情况。

截断奇异值正则化方法对直接对投影系数矩阵进行改造，寻找一个良态矩阵L_k，使其在二范数下较好的近似原始投影矩阵L。

$L_k=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{u}_i{\mathrm{\&phi;}}_i{v_i}^T---\left(10\right)$

使得

${\mathrm{\&Sigma;}}_{\mathrm{tsvd}}=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}\frac{<u_i,b>}{{\mathrm{\&phi;}}_i}{v}_i---\left(11\right)$

事实上，上式相当于将投影矩阵中数值较小的奇异值直接截断，从而有效的避免了数值较小的奇异值对散射噪声的放大。

全变分正则化在（6）式的基础上，加上了L1范数的J_TV(Σ)惩罚项，使得求解问题优化问题，

${\min \left|\right|\mathrm{L\&Sigma;}-B\left|\right|}_2^2+\mathrm{\&alpha;}{J}_{\mathrm{TV}}\left(\mathrm{\&Sigma;}\right)---\left(12\right)$

这里，

$J_{\mathrm{TV}}\left(\mathrm{\&Sigma;}\right)=\left|\right|\&dtri;\mathrm{\&Sigma;}\left|\right|={\&Integral;}_u\sqrt{{\mathrm{\&Sigma;}}_u^2+{\mathrm{\&beta;}}^2}\mathrm{du}---\left(13\right)$

表示全变差惩罚项，这里α指示正则参数，它表示了最小二乘项和惩罚项之间平衡比例。过大的正则参数和过小的正则化参数都不可取，一般通过一些数学准则确定正则参数的大小，这些准则包括广义交叉验证原理，Ｌ曲线方法等。

表示客体的线性吸收系数在水平方向的一阶导数。β是可调参数，避免（13）式中的不可微。

全变分的引入，对于重构客体相邻两种材料分界面有很好的刻画，与传统的Tikhonov方法相比，边缘保真度更高。

基于以上所述的截断奇异值和全变分正则化的方法以及所开发的闪光照相客体重建系统可精确实现被探测客体的线性吸收系数分布以及几何参数分布重建，并且具有健壮的抗噪声能力，满足工程要求。

Claims (6)

1.一种基于截断奇异值和全变分的闪光照相客体正则化重建方法，其特征在于，包括以下步骤：

（1）通过正向光子输运程序，获取闪光照相成像底片信息；

（2）通过截断奇异值方法获得闪光照相客体线性系数空间分布和几何参数；

（3）通过全变分正则化方法获得闪光照相客体线性系数空间分布和几何参数。

2.根据权利要求1所述的基于截断奇异值和全变分的闪光照相客体正则化重建方法，其特征在于，步骤1）具体包括以下步骤：

1）从成像底片上选取获得各个像素点的中心坐标；

2）光子输运模型计算得到成像底片像素点透射照射量信息；

3）从闪光成像底片透射照射量，利用闪光照相截断奇异值方法得到闪光照相客体吸收系数空间分布；

4）从闪光成像底片透射照射量，利用闪光照相全变差正则化方法得到闪光照相客体吸收系数空间分布。

3.根据权利要求2所述的基于截断奇异值和全变分的闪光照相客体正则化重建方法，其特征在于，所述的闪光照相截断奇异值方法为：从投影矩阵直接出发，获得投影矩阵的奇异值分解，并将投影矩阵数值小的奇异值截断，保留数值大的奇异值；实现对含散射噪声重构时误差放大的抑制。

4.根据权利要求2所述的基于截断奇异值和全变分的闪光照相客体正则化重建方法，其特征在于，所述光子输运模型为光子输运正向模型。

5.根据权利要求2所述的基于截断奇异值和全变分的闪光照相客体正则化重建方法，其特征在于，所述的闪光照相全变差正则化方法为：从投影的最小二乘方程出发，添加L1范数的惩罚项，并予以正则参数平衡最小二乘项和正则项的平衡关系；在获得高精度的重构结果的同时，具有边缘保真能力。

6.根据权利要求4所述的基于截断奇异值和全变分的闪光照相客体正则化重建方法，其特征在于，闪光照相正则参数的确立方法为：利用数学准则，最优化确定正则参数；所述数学准则包括广义交叉验证原理，Ｌ曲线原理，准最优化。

Images