CN113706411A - 一种高能闪光x射线图像非线性重建方法 - Google Patents

Abstract

本发明公开了图像处理技术领域的一种高能闪光X射线图像的重建方法，包括：获取X射线透射率图像；基于X射线透射率图像估算模糊核，构建非线性重建模型；基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；基于全条件概率密度函数，对后验概率密度函数进行MH抽样，获得样本值；采用截断牛顿共轭梯度法对样本值进行优化求解，获得MAP估计值。本发明能够有效提高X射线图像重建的精度，在医疗、工业等领域具有较高的工程应用价值和广阔的市场前景。

Description

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

技术领域

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

背景技术

X射线成像技术是研究核武器内部结构的重要手段。在X射线成像技术对高密度材料的诊断研究中，主要目标之一是准确测量客体内部空间密度分布。高能闪光X射线照相作为一项无损检测技术，可根据探测平面上X射线的空间强度分布实现照相目标空间密度分布的准确测量。但由于高能闪光X射线成像系统自身的复杂性，密度测量的精度容易受到系统模糊、散射以及噪声的影响。

马尔可夫链蒙特卡罗(Markov chain Monte Carlo,MCMC)方法作为一类随机方法，在求解高维反演问题中具有广泛的应用。目前，基于MCMC算法的高维反演算法主要通过构建线性重建模型进行求解。GPSR算法作为确定性方法可基于梯度投影技术对感兴趣的变量进行快速估计。LRIS(Gamma)和LRIS(Jeffreys) 重建算法分别基于Gamma先验和Jeffreys先验在线性模型下低秩近似Hessian 矩阵，并通过截断SVD计算目标参数的闭合解。线性重建忽略了系统模糊对于重建结果的影响，难以保证照相目标密度测量的精度。研究基于随机后优化和信赖域的X射线图像非线性重建方法，将开辟一条X射线图像重建的新途径，提高图像重建的精度，进而提高闪光实验中照相目标体密度分布测量的精度，对我国的国防建设也具有重要的研究价值和意义。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种高能闪光X射线图像非线性重建方法，能够解决系统模糊问题，提高X射线图像的重建精度。

为达到上述目的，本发明是采用下述技术方案实现的：

第一方面，本发明提供了一种高能闪光X射线图像非线性重建方法，包括：

获取X射线透射率图像；

基于X射线透射率图像估算模糊核，构建非线性重建模型；

基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；

基于全条件概率密度函数，对后验概率密度函数进行MH抽样，获得样本值；

采用截断牛顿共轭梯度法对样本值进行优化求解，获得MAP估计值。

进一步的，所述后验概率密度函数通过非线性重建模型引入超先验参数获得，包括：

所述非线性成像模型为：

y_t＝Kexp(-Hx)+n (21)

其中，y_t∈R^m为透射率图像的矢量形式，K∈R^m×m为系统模糊卷积核P_sys对应的矩阵形式，H∈R^m×n为成像系统的正向投影矩阵，n∈R^m为噪声项的矢量形式， n～N(0,ε^-1I)，其中ε为噪声精度参数；

采用Tikhonov正则化进行先验约束，将非线性重建模型转化成最小二乘问题为：

其中，

表示当括号内的函数值最小时线吸收系数x的取值，非线性算子F满足F(x)＝Kexp(-Hx)，L由GMRF(高斯马尔科夫随机场)定义；

表示F(x)与y的差值的二范数，其中F(x)＝K exp(-Hx)，y为透射率图像，

表示由GMRF定义的L与待求的线吸收系数x的差值的二范数，其中下标v是二范数

的参数，x_υ表示该式的估计值；

线吸收系数满足x～N(0,σ^-1L^ΤL)，其中超参数σ为先验精度参数；在贝叶斯框架下定义服从Gamma分布的超参数先验概率密度函数为p(σ)和p(ε)，则后验概率密度函数为：

其中，后验概率密度p(x,ε,σ|y)正比于以下四项的乘积：似然函数p(y|x,ε)，先验分布p(x|σ)，超参数先验概率密度函数为p(σ)和p(ε)；ε和σ为先验精度参数，α和β分别表示Gamma分布的形状参数与逆尺度参数，exp()表示自然对数e 的指数函数。

进一步的，所述全条件密度概率函数p(x|y,ε,σ)正比于以e为底数，以

为指数的指数函数，表示为：

先验精度参数ε和先验精度参数σ的条件密度均服从Gamma分布，表示为：

其中m和n分别表示成像系统的正向投影矩阵H∈R^m×n的行列数，α_ε和α_σ分别表示先验精度参数ε和σ对应的形状参数，β_ε和β_σ分别表示先验精度参数ε和σ对应的逆尺度参数。

进一步的，基于全条件概率密度函数，对后验概率密度函数进行MH抽样，包括：

定义增广正向模型的矩阵形式为：

其中，ε和σ为先验精度参数，F为非线性算子，满足F(x)＝K exp(-Hx)，L 由GMRF定义，x为待求的线吸收系数；

定义观测数据的矩阵形式为：

其中，ε为先验精度参数，y为透射率图像；

基于增广正向模型和观测数据计算最大后验MAP估计值，公式为：

其中，函数

表示二范数

取最小值时线吸收系数x的取值，ψ(x)为泛函；

基于增广正向模型和观测数据，求解以下随机优化问题得到RTO样本，即：

函数

表示二范

取最小值时泛函ψ的取值，其中

为稀疏QR分解J_ε,σ(x_ε,σ)＝Q_ε,σR_ε,σ中的Q_ε,σ的转置矩阵。

进一步的，所述重建方法通过信赖域方法将最大后验MAP估计值的求解转化为其子问题的求解并进行迭代，包括：

选定信赖域半径Δ_k＞0，通过二次近似模型构造信赖域子问题中的目标函数为：

式中，grad_k(x_k)和Hess_k(x_k)分别表示泛函

在第k次迭代点x_k处的梯度和Hessian矩阵，ξ为信赖域子问题中的自变量；

用F_ε,σ′(x_k)表示雅可比矩阵，计算grad_k(x_k)和Hess_k(x_k)为：

grad_k＝(F_ε,σ′(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ) (32)

Hess_k＝(F_ε,σ′(x_k))^Τ(F_ε,σ′(x_k))+(F_ε,σ″(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ) (33)

采用非精确的求解方法来求解泛函ψ的Hessian矩阵，忽略式(13)中第二项 (F_ε,σ″(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ)在真实解附近的值，以近似Hessian矩阵

代替原Hessian矩阵，信赖域子问题修改为：

其中，Hess_k(x_k)表示泛函在第k次迭代点x_k处的Hessian矩阵，用非精确求解方法求解，

表示近似Hessian矩阵；

在求解过程中，定义τ_k来表示目标模型的真实下降量Tard_k与近似模型的预估下降量Appd_k的比值，该判断准则下的信赖域方法要求每步迭代后的目标函数有所下降，此参数用于保证目标函数具有严格单调下降性，表示为：

其中Tard_k＝ψ(x_k)-ψ(x_k+ξ_k)，Appd_k＝φ_k(0)-φ_k(ξ_k)＝-φ_k(ξ_k)。

进一步的，采用截断牛顿共轭梯度法对样本值进行优化求解，包括：

结合截断牛顿共轭梯度法对上述信赖域子问题进行优化求解，牛顿共轭梯度法通过目标函数的梯度逐步生成共轭方向，并将获得的方向组作为线搜索方向实现二次函数最小值的求解，以初始方向作为负梯度方向，即：

d₀＝-grad₀ (36)

其中，-grad₀表示负梯度方向，d₀表示初始方向；

利用截断牛顿共轭梯度法进行信赖域子问题的非精确求解，生成的点列为：

ξ_l+1＝ξ_l+α_l+1d_l+1 (37)

式中，ξ为信赖域子问题中的自变量，α表示Gamma分布的形状参数，d为梯度方向，下标l为迭代次数；式中，

表示第l+1次梯度下降迭代目标函数的梯度，

表示第l+1次和第一次梯度下降迭代中，目标函数梯度的二范数的比值；

其中，第一次梯度下降迭代的目标函数梯度计算方法如下：

式中，

为Hessian矩阵近似值与第一次梯度下降迭代的ξ值的乘积， k表示迭代次数，grad_k为泛函第k次迭代的梯度；d_l和d_l ^Τ表示第一次梯度下降迭代的梯度方向及其转置矩阵，

表示Hessian矩阵近似值，

表示第一次梯度下降迭代的目标函数梯度的转置，求解过程中，满足

第二方面，一种高能闪光X射线图像非线性重建装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据上述任一项所述方法的步骤。

第三方面，计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述任一项所述方法的步骤。

与现有技术相比，本发明所达到的有益效果：

本发明在图像重建过程中将系统模糊考虑在内，解决了图像系统模糊对图像重建造成的困难，有效提升了图像重建的精度；本发明用基于信赖域技术的方法对模型进行优化求解，有效防止了过度求解，使得图像重建结果更佳；本发明采用截断牛顿共轭梯度法对信赖域子问题进行优化求解，有效提高了求解效率，降低计算成本。

附图说明

图1是本发明实施例一提供的高能闪光X射线图像非线性重建方法流程图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

实施例一：

一种X射线图像非线性重建方法，首先获取X射线透射率图像，结合由X 射线投影图像估计出的模糊核，这里采用基于稀疏先验L0正则化的模糊核估计，引入Tikhonov正则项，构建非线性重建模型；利用随机后优化的方法进行MH 抽样，并采用基于信赖域的非精确牛顿共轭梯度方法对抽样结果进行优化求解，具体步骤如下：

步骤1：假设图像存在噪声与系统模糊，构建透射率图像的非线性重建模型，并引入超先验参数，得到后验概率密度函数，具体步骤如下：

(11)假设图像存在噪声与系统模糊，非线性成像模型表示为：

y_t＝K exp(-Hx)+n (41)

其中，y_t∈R^m为透射率图像的矢量形式，K∈R^m×m为系统模糊卷积核P_sys对应的矩阵形式，H∈R^m×n为成像系统的正向投影矩阵，n∈R^m为噪声项的矢量形式， n～N(0,ε^-1I)，其中ε为噪声精度参数。

(12)采用Tikhonov正则化进行先验约束，将非线性重建模型转化成最小二乘问题为：

其中，

表示当括号内的函数值最小时线吸收系数x的取值，非线性算子F满足F(x)＝K exp(-Hx)，L由GMRF(高斯马尔科夫随机场)定义。

表示F(x)与y的差值的二范数，其中F(x)＝Kexp(-Hx)，y为透射率图像，

表示由GMRF(高斯马尔科夫随机场)定义的L与待求的线吸收系数x的差值的二范数，其中下标v是二范数

的参数，x_υ表示该式的估计值。

(13)假设线吸收系数满足x～N(0,σ^-1L^ΤL)，其中超参数σ为先验精度参数。在贝叶斯框架下定义服从Gamma分布的超参数先验概率密度函数为p(σ)和p(ε)，则后验概率密度函数为：

其中，后验概率密度p(x,ε,σ|y)正比于以下四项的乘积：似然函数p(y|x,ε)，先验分布p(x|σ)，超参数先验概率密度函数为p(σ)和p(ε)；ε和σ为先验精度参数，α和β分别表示Gamma分布的形状参数与逆尺度参数，exp()表示自然对数e 的指数函数，

和

的含义与式(2)相同。

步骤2：获得非线性重建模型的全条件概率密度函数，为了从后验概率中采样，需要得到全条件概率密度函数，二者公式中参数

相同，具体步骤如下：

(21)全条件密度概率函数p(x|y,ε,σ)正比于以e为底数，以

为指数的指数函数，表示为：

其中，

和

的含义与式(2)相同。

(22)先验精度参数ε和先验精度参数σ的条件密度均服从Gamma分布，表示为：

其中m和n分别表示式(1)中成像系统的正向投影矩阵H∈R^m×n的行列数，α_ε和α_σ分别表示先验精度参数ε和σ对应的形状参数，β_ε和β_σ分别表示先验精度参数ε和σ对应的逆尺度参数，二范数

和

的定义与式(2)相同。

步骤3：基于RTO(随机后优化)-MH(Metropolis-Hastings)算法，并采用基于信赖域技术的非精确牛顿共轭梯度法对样本值(即RTO样本)进行优化求解，即构建信赖域子问题的目标函数进行迭代求解，该步骤提供了一个优化的MAP 估计值求解方法，具体步骤如下：

(31)首先定义增广正向模型的矩阵形式为：

其中，ε和σ为先验精度参数，F为非线性算子，满足F(x)＝K exp(-Hx)，L 由GMRF(高斯马尔科夫随机场)定义，x为待求的线吸收系数。

(32)定义观测数据的矩阵形式为：

其中，ε为先验精度参数，y为透射率图像。

(33)在RTO方法的实现中，需要计算最大后验MAP估计值，基于公式(7) 和(8)计算MAP估计值，公式为：

其中，函数

表示二范数

取最小值时线吸收系数x的取值，其中F_ε,σ(x)和y_ε,σ由式(7)和式(8)定义。ψ(x)为泛函。

(34)在实际求解时，需要通过求解以下随机优化问题得到RTO样本，即 MH抽样：

函数

表示二范

取最小值时泛函ψ的取值，其中

为稀疏QR分解J_ε,σ(x_ε,σ)＝Q_ε,σR_ε,σ中的Q_ε,σ的转置矩阵，

和y_ε,σ由式(7)和式(8)定义。

(35)结合牛顿正则化方法局部二次快速收敛的特性、信赖域全局化对病态问题求解的有效性以及Eisenstat-Walker防止过度求解的思想，采用基于信赖域技术的非精确牛顿共轭梯度法(TRINCG)进行优化求解。

(36)信赖域方法将对式(9)中MAP的求解转化为其子问题的求解并进行迭代。首先寻找一个与式(9)相似的模型，选定信赖域半径Δ_k＞0，通过二次近似模型构造信赖域子问题中的目标函数为：

式中，grad_k(x_k)和Hess_k(x_k)分别表示泛函

在第k次迭代点x_k处的梯度和Hessian矩阵，ξ为信赖域子问题中的自变量。

用F_ε,σ′(x_k)表示雅可比矩阵，计算grad_k(x_k)和Hess_k(x_k)为：

grad_k＝(F_ε,σ′(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ) (52)

Hess_k＝(F_ε,σ′(x_k))^Τ(F_ε,σ′(x_k))+(F_ε,σ″(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ) (53)

构建子问题目标函数的目的是通过迭代，求解最大后验MAP估计值。

(37)采用非精确的求解方法来求解泛函ψ的Hessian矩阵，忽略式(13) 中第二项(F_ε,σ″(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ)在真实解附近的值，以近似Hessian矩阵

代替原Hessian矩阵，则信赖域子问题可以修改为：

其中，式(11)中的第二项中Hess_k(x_k)表示泛函在第k次迭代点x_k处的Hessian 矩阵，用非精确求解方法求解，所以式(14)中

的表示近似Hessian 矩阵。

在求解过程中，定义τ_k来表示目标模型的真实下降量Tard_k与近似模型的预估下降量Appd_k的比值，该判断准则下的信赖域方法要求每步迭代后的目标函数有所下降，此参数用于保证目标函数具有严格单调下降性，表示为：

其中Tard_k＝ψ(x_k)-ψ(x_k+ξ_k)，Appd_k＝φ_k(0)-φ_k(ξ_k)＝-φ_k(ξ_k)。

步骤4：采用截断牛顿共轭梯度法对信赖域子问题进行优化求解，具体步骤如下：

(41)结合截断牛顿共轭梯度法对上述信赖域子问题进行优化求解。牛顿共轭梯度法通过目标函数的梯度逐步生成共轭方向，并将获得的方向组作为线搜索方向实现二次函数最小值的求解。一般以初始方向作为负梯度方向，即：

d₀＝-grad₀ (56)

其中，-grad₀表示负梯度方向，d₀表示初始方向。

(42)利用截断牛顿共轭梯度法进行信赖域子问题(14)的非精确求解，生成的点列为：

ξ_l+1＝ξ_l+α_l+1d_l+1 (57)

式(17)中，ξ为信赖域子问题中的自变量，α表示Gamma分布的形状参数， d为梯度方向，下标l为迭代次数；式(18)中，

表示第l+1次梯度下降迭代目标函数的梯度，

表示第l+1次和第一次梯度下降迭代中，目标函数梯度的二范数的比值；

其中，

式(19)表示第一次梯度下降迭代的目标函数梯度计算方法，第一项

为Hessian矩阵近似值与第一次梯度下降迭代的ξ值的乘积，k表示迭代次数， grad_k为泛函第k次迭代的梯度(每次迭代中有梯度下降迭代，用l和l+1表示)；式(20)中d_l和d_l ^Τ表示第一次梯度下降迭代的梯度方向及其转置矩阵，

表示Hessian矩阵近似值，

表示第一次梯度下降迭代的目标函数梯度的转置。

求解过程中，不要求近似Hessian矩阵

的正定性，只需要满足

则可以顺利进行牛顿共轭梯度的迭代过程。

迭代获得的结果是MAP估计值，用于求解式(10)中的x^*，因为该求解方法是建立在含有模糊核的非线性重建模型基础上的，所以求解得到的MAP估计值以及进一步得到的线吸收系数RTO样本值x^*都有更高的精度，线吸收系数即为图像重建的结果，所以提高了线吸收系数的求解精度即为提高了图像重建精度。

实施例二：

本发明实施例还提供了一种高能闪光X射线图像非线性重建装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行下述方法的步骤：

获取X射线透射率图像；

基于X射线透射率图像估算模糊核，构建非线性重建模型；

基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；

基于全条件概率密度函数，对后验概率密度函数进行MH抽样，获得样本值；

采用截断牛顿共轭梯度法对样本值进行优化求解，获得MAP估计值。

实施例三：

本发明实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现下述方法的步骤：

获取X射线透射率图像；

基于X射线透射率图像估算模糊核，构建非线性重建模型；

基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；

基于全条件概率密度函数，对后验概率密度函数进行MH抽样，获得样本值；

采用截断牛顿共轭梯度法对样本值进行优化求解，获得MAP估计值。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、 CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/ 或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (8)

1.一种高能闪光X射线图像非线性重建方法，其特征是，包括：

获取X射线透射率图像；

基于X射线透射率图像估算模糊核，构建非线性重建模型；

基于非线性重建模型，得到后验概率密度函数和全条件概率密度函数；

基于全条件概率密度函数，对后验概率密度函数进行MH抽样，获得样本值；

采用截断牛顿共轭梯度法对样本值进行优化求解，获得MAP估计值。

2.根据权利要求1所述的高能闪光X射线图像非线性重建方法，其特征是，所述后验概率密度函数通过非线性重建模型引入超先验参数获得，包括：

所述非线性成像模型为：

y_t＝Kexp(-Hx)+n (1)

其中，y_t∈R^m为透射率图像的矢量形式，K∈R^m×m为系统模糊卷积核P_sys对应的矩阵形式，H∈R^m×n为成像系统的正向投影矩阵，n∈R^m为噪声项的矢量形式，n～N(0,ε^-1I)，其中ε为噪声精度参数；

采用Tikhonov正则化进行先验约束，将非线性重建模型转化成最小二乘问题为：

其中，

表示当括号内的函数值最小时线吸收系数x的取值，非线性算子F满足F(x)＝Kexp(-Hx)，L由GMRF(高斯马尔科夫随机场)定义；

表示F(x)与y的差值的二范数，其中F(x)＝Kexp(-Hx)，y为透射率图像，

表示由GMRF定义的L与待求的线吸收系数x的差值的二范数，其中下标v是二范数

的参数，x_υ表示该式的估计值；

线吸收系数满足x～N(0,σ^-1L^ΤL)，其中超参数σ为先验精度参数；在贝叶斯框架下定义服从Gamma分布的超参数先验概率密度函数为p(σ)和p(ε)，则后验概率密度函数为：

其中，后验概率密度p(x,ε,σ|y)正比于以下四项的乘积：似然函数p(y|x,ε)，先验分布p(x|σ)，超参数先验概率密度函数为p(σ)和p(ε)；ε和σ为先验精度参数，α和β分别表示Gamma分布的形状参数与逆尺度参数，exp()表示自然对数e的指数函数。

3.根据权利要求2所述的高能闪光X射线图像非线性重建方法，其特征是，所述全条件密度概率函数p(x|y,ε,σ)正比于以e为底数，以

为指数的指数函数，表示为：

先验精度参数ε和先验精度参数σ的条件密度均服从Gamma分布，表示为：

其中m和n分别表示成像系统的正向投影矩阵H∈R^m×n的行列数，α_ε和α_σ分别表示先验精度参数ε和σ对应的形状参数，β_ε和β_σ分别表示先验精度参数ε和σ对应的逆尺度参数。

4.根据权利要求3所述的高能闪光X射线图像非线性重建方法，其特征是，基于全条件概率密度函数，对后验概率密度函数进行MH抽样，包括：

定义增广正向模型的矩阵形式为：

其中，ε和σ为先验精度参数，F为非线性算子，满足F(x)＝Kexp(-Hx)，L由GMRF定义，x为待求的线吸收系数；

定义观测数据的矩阵形式为：

其中，ε为先验精度参数，y为透射率图像；

基于增广正向模型和观测数据计算最大后验MAP估计值，公式为：

其中，函数

表示二范数

取最小值时线吸收系数x的取值，ψ(x)为泛函；

基于增广正向模型和观测数据，求解以下随机优化问题得到RTO样本，即：

函数

表示二范

取最小值时泛函ψ的取值，其中

为稀疏QR分解J_ε,σ(x_ε,σ)＝Q_ε,σR_ε,σ中的Q_ε,σ的转置矩阵。

5.根据权利要求4所述的高能闪光X射线图像非线性重建方法，其特征是，所述重建方法通过信赖域方法将最大后验MAP估计值的求解转化为其子问题的求解并进行迭代，包括：

选定信赖域半径Δ_k＞0，通过二次近似模型构造信赖域子问题中的目标函数为：

式中，grad_k(x_k)和Hess_k(x_k)分别表示泛函

在第k次迭代点x_k处的梯度和Hessian矩阵，ξ为信赖域子问题中的自变量；

用F_ε,σ′(x_k)表示雅可比矩阵，计算grad_k(x_k)和Hess_k(x_k)为：

grad_k＝(F_ε,σ′(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ) (12)

Hess_k＝(F_ε,σ′(x_k))^Τ(F_ε,σ′(x_k))+(F_ε,σ″(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ) (13)

采用非精确的求解方法来求解泛函ψ的Hessian矩阵，忽略式(13)中第二项(F_ε,σ″(x_k))^Τ(F_ε,σ(x_k)-y_ε,σ)在真实解附近的值，以近似Hessian矩阵

代替原Hessian矩阵，信赖域子问题修改为：

其中，Hess_k(x_k)表示泛函在第k次迭代点x_k处的Hessian矩阵，用非精确求解方法求解，

表示近似Hessian矩阵；

在求解过程中，定义τ_k来表示目标模型的真实下降量Tard_k与近似模型的预估下降量Appd_k的比值，该判断准则下的信赖域方法要求每步迭代后的目标函数有所下降，此参数用于保证目标函数具有严格单调下降性，表示为：

其中Tard_k＝ψ(x_k)-ψ(x_k+ξ_k)，Appd_k＝φ_k(0)-φ_k(ξ_k)＝-φ_k(ξ_k)。

6.根据权利要求5所述的高能闪光X射线图像非线性重建方法，其特征是，采用截断牛顿共轭梯度法对样本值进行优化求解，包括：

结合截断牛顿共轭梯度法对上述信赖域子问题进行优化求解，牛顿共轭梯度法通过目标函数的梯度逐步生成共轭方向，并将获得的方向组作为线搜索方向实现二次函数最小值的求解，以初始方向作为负梯度方向，即：

d₀＝-grad₀ (16)

其中，-grad₀表示负梯度方向，d₀表示初始方向；

利用截断牛顿共轭梯度法进行信赖域子问题的非精确求解，生成的点列为：

ξ_l+1＝ξ_l+α_l+1d_l+1 (17)

式中，ξ为信赖域子问题中的自变量，α表示Gamma分布的形状参数，d为梯度方向，下标l为迭代次数；式中，

表示第l+1次梯度下降迭代目标函数的梯度，

表示第l+1次和第一次梯度下降迭代中，目标函数梯度的二范数的比值；

其中，第一次梯度下降迭代的目标函数梯度计算方法如下：

式中，

为Hessian矩阵近似值与第一次梯度下降迭代的ξ值的乘积，k表示迭代次数，grad_k为泛函第k次迭代的梯度；d_l和d_l ^Τ表示第一次梯度下降迭代的梯度方向及其转置矩阵，

表示Hessian矩阵近似值，

表示第一次梯度下降迭代的目标函数梯度的转置，求解过程中，满足

7.一种高能闪光X射线图像非线性重建装置，其特征在于，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据权利要求1～6任一项所述方法的步骤。

8.计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现权利要求1～6任一项所述方法的步骤。

Images