一种高能闪光X射线图像非线性重建方法

一种高能闪光x射线图像非线性重建方法
技术领域
1.本发明涉及一种高能闪光x射线图像非线性重建方法，属于图像处理技术领域。


背景技术：

2.x射线成像技术是研究核武器内部结构的重要手段。在x射线成像技术对高密度材料的诊断研究中，主要目标之一是准确测量客体内部空间密度分布。高能闪光x射线照相作为一项无损检测技术，可根据探测平面上x射线的空间强度分布实现照相目标空间密度分布的准确测量。但由于高能闪光x射线成像系统自身的复杂性，密度测量的精度容易受到系统模糊、散射以及噪声的影响。
3.马尔可夫链蒙特卡罗(markov chain monte carlo,mcmc)方法作为一类随机方法，在求解高维反演问题中具有广泛的应用。目前，基于mcmc算法的高维反演算法主要通过构建线性重建模型进行求解。gpsr算法作为确定性方法可基于梯度投影技术对感兴趣的变量进行快速估计。lris(gamma)和lris(jeffreys) 重建算法分别基于gamma先验和jeffreys先验在线性模型下低秩近似hessian 矩阵，并通过截断svd计算目标参数的闭合解。线性重建忽略了系统模糊对于重建结果的影响，难以保证照相目标密度测量的精度。研究基于随机后优化和信赖域的x射线图像非线性重建方法，将开辟一条x射线图像重建的新途径，提高图像重建的精度，进而提高闪光实验中照相目标体密度分布测量的精度，对我国的国防建设也具有重要的研究价值和意义。


技术实现要素：

4.本发明的目的在于克服现有技术中的不足，提供一种高能闪光x射线图像非线性重建方法，能够解决系统模糊问题，提高x射线图像的重建精度。
5.为达到上述目的，本发明是采用下述技术方案实现的：
6.第一方面，本发明提供了一种高能闪光x射线图像非线性重建方法，包括：
7.获取x射线透射率图像；
8.基于x射线透射率图像估算模糊核，构建非线性重建模型；
9.基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；
10.基于全条件概率密度函数，对后验概率密度函数进行mh抽样，获得样本值；
11.采用截断牛顿共轭梯度法对样本值进行优化求解，获得map估计值。
12.进一步的，所述后验概率密度函数通过非线性重建模型引入超先验参数获得，包括：
13.所述非线性成像模型为：
14.y
t
＝kexp(
‑
hx)+n
ꢀꢀꢀ
(21)
15.其中，y
t
∈r
m
为透射率图像的矢量形式，k∈r
m
×
m
为系统模糊卷积核p
sys
对应的矩阵形式，h∈r
m
×
n
为成像系统的正向投影矩阵，n∈r
m
为噪声项的矢量形式， n～n(0,ε
‑1i)，其中ε为噪声精度参数；
16.采用tikhonov正则化进行先验约束，将非线性重建模型转化成最小二乘问题为：
[0017][0018]
其中，表示当括号内的函数值最小时线吸收系数x的取值，非线性算子f满足f(x)＝kexp(
‑
hx)，l由gmrf(高斯马尔科夫随机场)定义；表示f(x)与y的差值的二范数，其中f(x)＝k exp(
‑
hx)，y为透射率图像，表示由gmrf定义的l与待求的线吸收系数x的差值的二范数，其中下标v是二范数的参数，x
υ
表示该式的估计值；
[0019]
线吸收系数满足x～n(0,σ
‑1l
τ
l)，其中超参数σ为先验精度参数；在贝叶斯框架下定义服从gamma分布的超参数先验概率密度函数为p(σ)和p(ε)，则后验概率密度函数为：
[0020][0021]
其中，后验概率密度p(x,ε,σ|y)正比于以下四项的乘积：似然函数p(y|x,ε)，先验分布p(x|σ)，超参数先验概率密度函数为p(σ)和p(ε)；ε和σ为先验精度参数，α和β分别表示gamma分布的形状参数与逆尺度参数，exp()表示自然对数e 的指数函数。
[0022]
进一步的，所述全条件密度概率函数p(x|y,ε,σ)正比于以e为底数，以为指数的指数函数，表示为：
[0023][0024]
先验精度参数ε和先验精度参数σ的条件密度均服从gamma分布，表示为：
[0025][0026][0027]
其中m和n分别表示成像系统的正向投影矩阵h∈r
m
×
n
的行列数，α
ε
和α
σ
分别表示先验精度参数ε和σ对应的形状参数，β
ε
和β
σ
分别表示先验精度参数ε和σ对应的逆尺度参数。
[0028]
进一步的，基于全条件概率密度函数，对后验概率密度函数进行mh抽样，包括：
[0029]
定义增广正向模型的矩阵形式为：
[0030][0031]
其中，ε和σ为先验精度参数，f为非线性算子，满足f(x)＝k exp(
‑
hx)，l 由gmrf定义，x为待求的线吸收系数；
[0032]
定义观测数据的矩阵形式为：
[0033][0034]
其中，ε为先验精度参数，y为透射率图像；
[0035]
基于增广正向模型和观测数据计算最大后验map估计值，公式为：
[0036][0037]
其中，函数表示二范数取最小值时线吸收系数x的取值，ψ(x)为泛函；
[0038]
基于增广正向模型和观测数据，求解以下随机优化问题得到rto样本，即：
[0039][0040]
函数表示二范取最小值时泛函ψ的取值，其中为稀疏qr分解j
ε,σ
(x
ε,σ
)＝q
ε,σ
r
ε,σ
中的q
ε,σ
的转置矩阵。
[0041]
进一步的，所述重建方法通过信赖域方法将最大后验map估计值的求解转化为其子问题的求解并进行迭代，包括：
[0042]
选定信赖域半径δ
k
＞0，通过二次近似模型构造信赖域子问题中的目标函数为：
[0043][0044]
式中，grad
k
(x
k
)和hess
k
(x
k
)分别表示泛函在第k次迭代点x
k
处的梯度和hessian矩阵，ξ为信赖域子问题中的自变量；
[0045]
用f
ε,σ
′
(x
k
)表示雅可比矩阵，计算grad
k
(x
k
)和hess
k
(x
k
)为：
[0046]
grad
k
＝(f
ε,σ
′
(x
k
))
τ
(f
ε,σ
(x
k
)
‑
y
ε,σ
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(32)
[0047]
hess
k
＝(f
ε,σ
′
(x
k
))
τ
(f
ε,σ
′
(x
k
))+(f
ε,σ
″
(x
k
))
τ
(f
ε,σ
(x
k
)
‑
y
ε,σ
)
ꢀꢀꢀꢀꢀ
(33)
[0048]
采用非精确的求解方法来求解泛函ψ的hessian矩阵，忽略式(13)中第二项 (f
ε,σ
″
(x
k
))
τ
(f
ε,σ
(x
k
)
‑
y
ε,σ
)在真实解附近的值，以近似hessian矩阵代替原hessian矩阵，信赖域子问题修改为：
[0049][0050]
其中，hess
k
(x
k
)表示泛函在第k次迭代点x
k
处的hessian矩阵，用非精确求解方法求解，表示近似hessian矩阵；
[0051]
在求解过程中，定义τ
k
来表示目标模型的真实下降量tard
k
与近似模型的预估下降量appd
k
的比值，该判断准则下的信赖域方法要求每步迭代后的目标函数有所下降，此参数用于保证目标函数具有严格单调下降性，表示为：
[0052][0053]
其中tard
k
＝ψ(x
k
)
‑
ψ(x
k
+ξ
k
)，appd
k
＝φ
k
(0)
‑
φ
k
(ξ
k
)＝
‑
φ
k
(ξ
k
)。
[0054]
进一步的，采用截断牛顿共轭梯度法对样本值进行优化求解，包括：
[0055]
结合截断牛顿共轭梯度法对上述信赖域子问题进行优化求解，牛顿共轭梯度法通过目标函数的梯度逐步生成共轭方向，并将获得的方向组作为线搜索方向实现二次函数最小值的求解，以初始方向作为负梯度方向，即：
[0056]
d0＝
‑
grad0ꢀꢀꢀ
(36)
[0057]
其中，
‑
grad0表示负梯度方向，d0表示初始方向；
[0058]
利用截断牛顿共轭梯度法进行信赖域子问题的非精确求解，生成的点列为：
[0059]
ξ
l+1
＝ξ
l
+α
l+1
d
l+1
ꢀꢀꢀ
(37)
[0060][0061]
式中，ξ为信赖域子问题中的自变量，α表示gamma分布的形状参数，d为梯度方向，下标l为迭代次数；式中，表示第l+1次梯度下降迭代目标函数的梯度，表示第l+1次和第一次梯度下降迭代中，目标函数梯度的二范数的比值；
[0062]
其中，第一次梯度下降迭代的目标函数梯度计算方法如下：
[0063][0064][0065]
式中，为hessian矩阵近似值与第一次梯度下降迭代的ξ值的乘积， k表示迭代次数，grad
k
为泛函第k次迭代的梯度；d
l
和d
lτ
表示第一次梯度下降迭代的梯度方向及其转置矩阵，表示hessian矩阵近似值，表示第一次梯度下降迭代的目标函数梯度的转置，求解过程中，满足
[0066]
第二方面，一种高能闪光x射线图像非线性重建装置，包括处理器及存储介质；
[0067]
所述存储介质用于存储指令；
[0068]
所述处理器用于根据所述指令进行操作以执行根据上述任一项所述方法的步骤。
[0069]
第三方面，计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述任一项所述方法的步骤。
[0070]
与现有技术相比，本发明所达到的有益效果：
[0071]
本发明在图像重建过程中将系统模糊考虑在内，解决了图像系统模糊对图像重建造成的困难，有效提升了图像重建的精度；本发明用基于信赖域技术的方法对模型进行优化求解，有效防止了过度求解，使得图像重建结果更佳；本发明采用截断牛顿共轭梯度法对
信赖域子问题进行优化求解，有效提高了求解效率，降低计算成本。
附图说明
[0072]
图1是本发明实施例一提供的高能闪光x射线图像非线性重建方法流程图。
具体实施方式
[0073]
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
[0074]
实施例一：
[0075]
一种x射线图像非线性重建方法，首先获取x射线透射率图像，结合由x 射线投影图像估计出的模糊核，这里采用基于稀疏先验l0正则化的模糊核估计，引入tikhonov正则项，构建非线性重建模型；利用随机后优化的方法进行mh 抽样，并采用基于信赖域的非精确牛顿共轭梯度方法对抽样结果进行优化求解，具体步骤如下：
[0076]
步骤1：假设图像存在噪声与系统模糊，构建透射率图像的非线性重建模型，并引入超先验参数，得到后验概率密度函数，具体步骤如下：
[0077]
(11)假设图像存在噪声与系统模糊，非线性成像模型表示为：
[0078]
y
t
＝k exp(
‑
hx)+n
ꢀꢀꢀ
(41)
[0079]
其中，y
t
∈r
m
为透射率图像的矢量形式，k∈r
m
×
m
为系统模糊卷积核p
sys
对应的矩阵形式，h∈r
m
×
n
为成像系统的正向投影矩阵，n∈r
m
为噪声项的矢量形式， n～n(0,ε
‑1i)，其中ε为噪声精度参数。
[0080]
(12)采用tikhonov正则化进行先验约束，将非线性重建模型转化成最小二乘问题为：
[0081][0082]
其中，表示当括号内的函数值最小时线吸收系数x的取值，非线性算子f满足f(x)＝k exp(
‑
hx)，l由gmrf(高斯马尔科夫随机场)定义。表示f(x)与y的差值的二范数，其中f(x)＝kexp(
‑
hx)，y为透射率图像，表示由gmrf(高斯马尔科夫随机场)定义的l与待求的线吸收系数x的差值的二范数，其中下标v是二范数的参数，x
υ
表示该式的估计值。
[0083]
(13)假设线吸收系数满足x～n(0,σ
‑1l
τ
l)，其中超参数σ为先验精度参数。在贝叶斯框架下定义服从gamma分布的超参数先验概率密度函数为p(σ)和p(ε)，则后验概率密度函数为：
[0084][0085]
其中，后验概率密度p(x,ε,σ|y)正比于以下四项的乘积：似然函数p(y|x,ε)，先验
分布p(x|σ)，超参数先验概率密度函数为p(σ)和p(ε)；ε和σ为先验精度参数，α和β分别表示gamma分布的形状参数与逆尺度参数，exp()表示自然对数e 的指数函数，和的含义与式(2)相同。
[0086]
步骤2：获得非线性重建模型的全条件概率密度函数，为了从后验概率中采样，需要得到全条件概率密度函数，二者公式中参数相同，具体步骤如下：
[0087]
(21)全条件密度概率函数p(x|y,ε,σ)正比于以e为底数，以为指数的指数函数，表示为：
[0088][0089]
其中，和的含义与式(2)相同。
[0090]
(22)先验精度参数ε和先验精度参数σ的条件密度均服从gamma分布，表示为：
[0091][0092][0093]
其中m和n分别表示式(1)中成像系统的正向投影矩阵h∈r
m
×
n
的行列数，α
ε
和α
σ
分别表示先验精度参数ε和σ对应的形状参数，β
ε
和β
σ
分别表示先验精度参数ε和σ对应的逆尺度参数，二范数和的定义与式(2)相同。
[0094]
步骤3：基于rto(随机后优化)
‑
mh(metropolis
‑
hastings)算法，并采用基于信赖域技术的非精确牛顿共轭梯度法对样本值(即rto样本)进行优化求解，即构建信赖域子问题的目标函数进行迭代求解，该步骤提供了一个优化的map 估计值求解方法，具体步骤如下：
[0095]
(31)首先定义增广正向模型的矩阵形式为：
[0096][0097]
其中，ε和σ为先验精度参数，f为非线性算子，满足f(x)＝k exp(
‑
hx)，l 由gmrf(高斯马尔科夫随机场)定义，x为待求的线吸收系数。
[0098]
(32)定义观测数据的矩阵形式为：
[0099][0100]
其中，ε为先验精度参数，y为透射率图像。
[0101]
(33)在rto方法的实现中，需要计算最大后验map估计值，基于公式(7) 和(8)计算
map估计值，公式为：
[0102][0103]
其中，函数表示二范数取最小值时线吸收系数x的取值，其中f
ε,σ
(x)和y
ε,σ
由式(7)和式(8)定义。ψ(x)为泛函。
[0104]
(34)在实际求解时，需要通过求解以下随机优化问题得到rto样本，即 mh抽样：
[0105][0106]
函数表示二范取最小值时泛函ψ的取值，其中为稀疏qr分解j
ε,σ
(x
ε,σ
)＝q
ε,σ
r
ε,σ
中的q
ε,σ
的转置矩阵，和y
ε,σ
由式(7)和式(8)定义。
[0107]
(35)结合牛顿正则化方法局部二次快速收敛的特性、信赖域全局化对病态问题求解的有效性以及eisenstat
‑
walker防止过度求解的思想，采用基于信赖域技术的非精确牛顿共轭梯度法(trincg)进行优化求解。
[0108]
(36)信赖域方法将对式(9)中map的求解转化为其子问题的求解并进行迭代。首先寻找一个与式(9)相似的模型，选定信赖域半径δ
k
＞0，通过二次近似模型构造信赖域子问题中的目标函数为：
[0109][0110]
式中，grad
k
(x
k
)和hess
k
(x
k
)分别表示泛函在第k次迭代点x
k
处的梯度和hessian矩阵，ξ为信赖域子问题中的自变量。
[0111]
用f
ε,σ
′
(x
k
)表示雅可比矩阵，计算grad
k
(x
k
)和hess
k
(x
k
)为：
[0112]
grad
k
＝(f
ε,σ
′
(x
k
))
τ
(f
ε,σ
(x
k
)
‑
y
ε,σ
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(52)
[0113]
hess
k
＝(f
ε,σ
′
(x
k
))
τ
(f
ε,σ
′
(x
k
))+(f
ε,σ
″
(x
k
))
τ
(f
ε,σ
(x
k
)
‑
y
ε,σ
)
ꢀꢀꢀꢀꢀ
(53)
[0114]
构建子问题目标函数的目的是通过迭代，求解最大后验map估计值。
[0115]
(37)采用非精确的求解方法来求解泛函ψ的hessian矩阵，忽略式(13) 中第二项(f
ε,σ
″
(x
k
))
τ
(f
ε,σ
(x
k
)
‑
y
ε,σ
)在真实解附近的值，以近似hessian矩阵代替原hessian矩阵，则信赖域子问题可以修改为：
[0116][0117]
其中，式(11)中的第二项中hess
k
(x
k
)表示泛函在第k次迭代点x
k
处的hessian 矩阵，用非精确求解方法求解，所以式(14)中的表示近似hessian 矩阵。
[0118]
在求解过程中，定义τ
k
来表示目标模型的真实下降量tard
k
与近似模型的预估下降量appd
k
的比值，该判断准则下的信赖域方法要求每步迭代后的目标函数有所下降，此参数用于保证目标函数具有严格单调下降性，表示为：
[0119][0120]
其中tard
k
＝ψ(x
k
)
‑
ψ(x
k
+ξ
k
)，appd
k
＝φ
k
(0)
‑
φ
k
(ξ
k
)＝
‑
φ
k
(ξ
k
)。
[0121]
步骤4：采用截断牛顿共轭梯度法对信赖域子问题进行优化求解，具体步骤如下：
[0122]
(41)结合截断牛顿共轭梯度法对上述信赖域子问题进行优化求解。牛顿共轭梯度法通过目标函数的梯度逐步生成共轭方向，并将获得的方向组作为线搜索方向实现二次函数最小值的求解。一般以初始方向作为负梯度方向，即：
[0123]
d0＝
‑
grad0ꢀꢀꢀ
(56)
[0124]
其中，
‑
grad0表示负梯度方向，d0表示初始方向。
[0125]
(42)利用截断牛顿共轭梯度法进行信赖域子问题(14)的非精确求解，生成的点列为：
[0126]
ξ
l+1
＝ξ
l
+α
l+1
d
l+1
ꢀꢀꢀ
(57)
[0127][0128]
式(17)中，ξ为信赖域子问题中的自变量，α表示gamma分布的形状参数， d为梯度方向，下标l为迭代次数；式(18)中，表示第l+1次梯度下降迭代目标函数的梯度，表示第l+1次和第一次梯度下降迭代中，目标函数梯度的二范数的比值；
[0129]
其中，
[0130][0131][0132]
式(19)表示第一次梯度下降迭代的目标函数梯度计算方法，第一项为hessian矩阵近似值与第一次梯度下降迭代的ξ值的乘积，k表示迭代次数， grad
k
为泛函第k次迭代的梯度(每次迭代中有梯度下降迭代，用l和l+1表示)；式(20)中d
l
和d
lτ
表示第一次梯度下降迭代的梯度方向及其转置矩阵，表示hessian矩阵近似值，表示第一次梯度下降迭代的目标函数梯度的转置。
[0133]
求解过程中，不要求近似hessian矩阵的正定性，只需要满足则可以顺利进行牛顿共轭梯度的迭代过程。
[0134]
迭代获得的结果是map估计值，用于求解式(10)中的x
*
，因为该求解方法是建立在含有模糊核的非线性重建模型基础上的，所以求解得到的map估计值以及进一步得到的线吸收系数rto样本值x
*
都有更高的精度，线吸收系数即为图像重建的结果，所以提高了线吸收系数的求解精度即为提高了图像重建精度。
[0135]
实施例二：
[0136]
本发明实施例还提供了一种高能闪光x射线图像非线性重建装置，包括处理器及存储介质；
[0137]
所述存储介质用于存储指令；
[0138]
所述处理器用于根据所述指令进行操作以执行下述方法的步骤：
[0139]
获取x射线透射率图像；
[0140]
基于x射线透射率图像估算模糊核，构建非线性重建模型；
[0141]
基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；
[0142]
基于全条件概率密度函数，对后验概率密度函数进行mh抽样，获得样本值；
[0143]
采用截断牛顿共轭梯度法对样本值进行优化求解，获得map估计值。
[0144]
实施例三：
[0145]
本发明实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现下述方法的步骤：
[0146]
获取x射线透射率图像；
[0147]
基于x射线透射率图像估算模糊核，构建非线性重建模型；
[0148]
基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；
[0149]
基于全条件概率密度函数，对后验概率密度函数进行mh抽样，获得样本值；
[0150]
采用截断牛顿共轭梯度法对样本值进行优化求解，获得map估计值。
[0151]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、 cd
‑
rom、光学存储器等)上实施的计算机程序产品的形式。
[0152]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/ 或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0153]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0154]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0155]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人
员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。