CN112053307A - 一种x射线图像线性重建方法 - Google Patents

Abstract

本发明公开了一种X射线图像线性重建方法，首先获取X射线投影图像，引入基于全变差先验的正则项，构造重建目标函数；引入超先验参数，构造分层贝叶斯模型；运用变量分裂方法，引入分裂变量，分离数据保真项和TV正则项，得到分裂形式下的联合概率密度分布；基于Jefferys先验定义超先验变量，得到各变量的条件分布；迭代更新超先验参数，对含TV正则项的分裂变量的条件分布进行求解；利用正向矩阵的低秩性质，近似待求参数的全条件概率密度分布，计算低秩近似的目标分布得到关于待求参数的闭合解；计算采样样本的均值，对待求参数进行估计。本发明能够有效地解决求解大规模线性逆问题中存在的计算开销大等问题。

Description

一种X射线图像线性重建方法

技术领域

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

背景技术

X射线成像技术是研究核武器内部结构的重要手段，是获取非核内爆演化过程后期弹芯物理特性和几何特性的主要工具。在X射线成像技术对高密度材料的诊断研究中，主要目标有两个：一是准确测量客体内部空间密度及其分布，二是定量确定客体内部的几何界面。客体的密度和界面测量均属于典型的高维反演问题，而且存在数据维数高，易受噪声、散射、光源及探测器模糊等系统模糊影响的难题。图像重建技术是X射线图像处理中的重要研究内容，先前针对图像重建的方法有滤波反投影重建法(FBP)、代数法(ART)和约束共轭梯度重建法(CCG)等，但存在重建精度低且对噪声的抑制能力较差等问题。

为了缓解重建反演问题的病态性，从低信噪比的投影图像中获得高精度的客体信息，近年来发展了多种图像重建方法，包括变分优化方法以及基于马尔可夫链蒙特卡罗(MCMC)的随机采样方法。变分优化方法在计算效率方面具有显著的优势，如梯度投影稀疏重建(GPSR)方法和改进的GPSR-BB方法，此类方法能对感兴趣的变量进行快速的估计，但不能提供参数的不确定度分析，且需要手动调整正则化参数。MCMC方法作为一类随机方法，在求解高维反演问题中具有广泛的应用。它能够对需要估计的参数进行不确定性量化，其存在的瓶颈在于高维高斯随机变量的重复采样，且在每次迭代时都需要进行矩阵的因式分解，计算成本较大。另一个挑战是，后验分布通常没有一个闭合解，在这种情况下，近似技术成为必要。同时，建立分层模型并引入变量分裂方法，将对提高X射线图像重建的高效性和精度有重要作用。研究基于分层模型和低秩近似的X射线图像重建方法，将开辟一条X射线图像重建的新途径，在X射线图像中的应用中具有明显优势。利用这种方法重建得到高质量的X射线图像，将提高闪光实验中客体界面和密度测量的精度，对我国的国防建设也具有重要的研究价值和意义。

发明内容

本发明所要解决的技术问题是：X射线图像重建计算成本易受噪声等因素影响，提供一种X射线图像线性重建方法，能够加速X射线图像重建进程，提高X射线图像的密度重建精度。

为解决上述技术问题，本发明提出了一种X射线图像线性重建方法，包括以下步骤：

步骤1：获取X射线投影图像y，引入基于全变差(TV)先验的正则项，构造重建的目标函数；

步骤2：引入超先验参数，构造分层贝叶斯模型，建立最大后验概率模型；所述超先验参数为在参数选择之前需要设置的参数，即将参数先验概率分布中的参数在贝叶斯推理过程中看作随机变量进行处理，可以看作是对于先验信息的一个预处理，所述参数先验概率分布中的参数包括均值、方差等；所述分层贝叶斯模型为先验概率密度函数含有某些超先验参数，将所述超先验参数视作随机变量，并服从某一先验分布，以提升参数选择的效果；

步骤3：运用变量分裂方法，引入分裂变量z，分离数据保真项和全变差TV正则项，得到分裂形式下的联合概率密度分布；所述变量分裂方法是引入分裂变量z，分离数据保真项和全变差TV正则项，然后求解相应的最小化问题；

步骤4：基于Jefferys先验定义超先验变量，得到各变量的条件分布；迭代更新超先验参数，对含全变差TV正则项的分裂变量z的条件分布进行求解；所述Jefferys先验，即非信息先验分布，其最主要性质就是不变性，即先验的形式不随着原分布参数形式的变化而变化；

步骤5：利用正向矩阵的低秩性质，近似变量x的全条件概率密度分布，计算低秩近似的目标分布得到关于变量x的闭合解，计算采样样本的均值，对变量x进行估计。

本发明所达到的有益效果：本发明提出了一种X射线图像线性重建方法，首先引入基于TV先验的正则项来构造重建的目标函数，能够在图像重建中有效地去除条形伪影并有效保留图像边缘信息，同时具有良好的抗噪性；引入超先验参数，构造分层贝叶斯模型，有利于提高参数获取的精度和效率；运用变量分裂方法，引入分裂变量z，分离数据保真项和TV正则项，得到分裂形式下的联合概率密度分布，简化了目标分布，有利于进行更简单的采样步骤，以便高效的先进的MCMC算法可以嵌入到每个采样任务；利用正向矩阵的低秩性质，近似图像x的全条件概率密度分布，以低秩近似作为采样的基础来加速高维高斯分布的绘制，能够有效地解决求解大规模线性逆问题时存在的计算开销大等问题。相较于其他方法，本发明可以加速X射线图像重建进程，更好地抑制测量数据和样本均值中存在的噪声和伪影，提高X射线图像的密度重建精度。

具体实施方式

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

本发明所描述的一种X射线图像线性重建方法，具体步骤如下：

步骤1：获取X射线投影图像y，引入基于全变差TV先验的正则项，构造重建的目标函数；

假设投影图像y被加性噪声破坏，则正向问题的线性随机模型为

y＝Hx+e (1)

其中，x∈R^N为需要重建的客体密度分布；y∈R^M为观测的X射线投影图像；H∈R^M×N为数值离散后得到的系统正向矩阵，又称为参数-观测数据映射矩阵；e∈R^M表示加性噪声或者测量误差，在本发明中，假设e是一个均值为0、方差为1/μ的高斯随机变量，与未知的变量x独立不相关，这里μ表示噪声参数，R表示实数空间，R^N和R^M×N分别表示N维和M×N维的实数空间，x∈R^N，y∈R^M说明x,y分别为N×1和M×1的列向量，H∈R^M×N说明H为M×N的矩阵。

TV正则项已经成为一种普遍存在的正则项，可用于解决大规模图像逆问题。

||x||_TV＝∑_1≤i,j≤N||(▽x)_i,j||₂ (2)

其中，||·||₂表示向量的二范数，即向量元素绝对值的平方和再开方，▽x表示变量x的二维离散梯度算子，在此设定下，关于变量x的后验概率密度分布p(x|y)为

其中，exp[·]表示以e为底数的指数函数，∝表示成比例，基于线性模型式(1)，从观测投影图像y∈R^M中重建出变量x∈R^N是一个病态逆问题，因此，采用TV正则化方法建模，并最小化如下的重建目标函数进行求解：

因此，重建的目标函数是对变量x的条件后验函数的最大化，也就是对函数

的最小化，其中，第一项为数据保真项；第二项为TV正则项；

为重建的客体密度分布；β为正则化参数，用来控制保真项与正则项之间的加权比例。

步骤2：引入超先验参数，构造分层贝叶斯模型，建立最大后验概率模型；所述超先验参数为在参数选择之前需要设置的参数，即将参数先验概率分布中的均值、方差等参数在贝叶斯推理过程中看作随机变量进行处理，可以看作是对于先验信息的一个预处理；所述分层贝叶斯模型，即先验概率密度函数中含有某些超参数，将所述超参数视作随机变量，并令随机变量服从某一先验分布，以提升参数选择的效果。

最大后验概率模型是根据经验数据对难以观察的量进行估计，应用贝叶斯理论建立最大后验概率模型，从定义的后验密度函数进行采样，去实现p(x|y)的最大化。首先，基于线性模型式(1)，在已知变量x,μ的条件下求投影图像y可表示为y|x,μ～N(Hx,μ^-1I)，将投影图像作为模型的一个变量，其中，～表示服从某一分布，N(·)表示高斯分布，I表示单位矩阵，于是可得到似然函数p(y|x,μ)为：

关于变量x的先验分布编码了在考虑观测数据之前期望或希望在x上执行的结构。假设变量x的先验分布服从高斯分布，

为了不事先选择最优参数，本发明中将噪声参数μ和正则化参数β作为未知数，视为随机变量，并给噪声参数μ和正则化参数β分配一个先验分布p(μ,β)，则噪声参数μ、正则化参数β和其他参数一起被估计，此时，噪声参数μ和正则化参数β亦被称为精度变量，在这种情况下，联合后验概率密度函数为：

由于，噪声参数μ和正则化参数β具有先验分布，联合后验概率密度函数不再服从高斯分布，不能以闭合解的形式存在，当p(μ)和p(β)为条件共轭先验时，变量x与噪声参数μ和正则化参数β分开进行迭代采样，通常情况下，噪声参数μ和正则化参数β单独更新。

步骤3：运用变量分裂方法，引入分裂变量z，分离数据保真项和TV正则项，得到分裂形式下的联合概率密度分布；

在贝叶斯推理框架内，MCMC方法的优点是对待求参数的后验分布提供了全面的描述，但在高维问题中，MCMC方法存在计算量大的瓶颈，为了克服这个限制，将变量分裂的方法嵌入到分层贝叶斯模型中，并推导出相应的具有封闭形式的条件后验分布，有利于解决大规模贝叶斯推理问题。

所述变量分裂方法是引入分裂变量z，分离数据保真项和全变差TV正则项，然后求解相应的最小化问题，

其中，f()代表数据保真项；r()代表一些正则化函数，arg表示对函数求参数(集合)，subject to表示约束条件。r()通常是非光滑的，甚至是非凸的；等式约束x＝z确保了求解式(8)等价于求解初始目标分布式(4)。

在贝叶斯推理问题中，变量分裂的目的是在一个优化子问题中单独使用目标函数的数据保真项f()和正则化函数r()，这种分而治之的策略有利于产生更简单的条件分布，从而更容易进行采样。引入一个正常数ρ来平衡变量x和分裂变量z之间的相似性，然后定义联合概率分布p(x,z|ρ)为

其中，

是一个散度函数，其作用是保证联合概率分布p(x,z|ρ)是合理的。在本发明中，散度函数

是二次的，并满足

根据式(9)可知，关于变量x和分裂变量z的条件概率分布为

结合式(7)，更新后的联合概率密度分布为

步骤4：基于Jefferys先验定义超先验变量，得到各变量的条件分布；迭代更新超先验参数，对含TV正则项的分裂变量z的条件分布进行求解；所述Jefferys先验，即非信息先验分布，其最主要性质就是不变性，即先验的形式不随着原分布参数形式的变化而变化。

本发明基于Jeffreys先验来定义噪声参数μ和正则化参数β，为了方便起见，利用方差变量而不是精度变量来进行模型参数化，即κ²＝μ^-1，τ²＝β^-1，于是，关于(κ²,τ²)的Jeffreys先验形式为

因此，正如Jeffreys所提倡的那样，κ²为尺度不变先验，τ²为数据水平方差尺度先验。在这种情况下，为了获得所需的全条件密度分布函数，定义精度变量ν＝τ²/κ²，进行参数变换后，步骤3中分裂形式下的联合后验概率密度函数(12)更新为

于是，可推导出各个变量的全条件概率密度函数为

其中，Γ(·)表示Gamma分布，因此，有

ν|x,z,κ²,y～1/Γ((N+2)/2,||z||_TV/(2κ²))，并可以依据以上分布迭代更新超先验参数。ν|x,z,κ²,y表示已知变量x,z,κ²,y的条件下求变量ν。

至于分裂变量z的求解，由于TV先验的非共轭性以及不可微性，可利用proximalMoreau-Yoshida-unadjusted Langevin算法(P-MYULA)(参照文献：Efficient Bayesiancomputation by proximal Markov Chain Monte Carlo:When Langevin meets Moreau)有效地对关于分裂变量z的全条件概率分布(18)进行高维采样。

步骤5：利用正向矩阵的低秩性质，近似图像x的全条件概率密度分布，计算低秩近似的目标分布得到关于变量x的闭合解，计算采样样本的均值，得到参数x的估计值。

在全贝叶斯框架下，后验分布通常没有一个闭合解，在这种情况下，可采用间接的基于采样的近似方法来研究后验概率密度分布，同时，在大规模的线性逆问题中，由于在每次迭代时都需要进行协方差矩阵的因式分解，同时高维高斯随机变量的重复采样也带来巨大的计算负担。为了克服这些问题，本发明利用正向矩阵的低秩结构对关于变量x的全条件密度分布进行低秩近似，使得目标分布便于采样。

由变量x的全条件密度分布(17)可得，

根据条件分布(19)进行采样得到的样本可以用x_s＝m_x+Gε来表示，其中随机变量ε～N(0,I)，ε为一个均值为0、方差为I高斯随机变量，条件协方差矩阵Γ_x的因式分解G满足Γ_x＝GG^T，计算均值m_x和随机向量Gε涉及协方差矩阵的因式分解，需要花费巨大的运算成本。利用正向矩阵H的低秩性质，可以构造一个可简易采样的快速目标分布，本发明中，通过一个截断的奇异值分解来近似H^TH，即

H^TH≈V_kΛ_kV_k ^T (20)

其中，酉矩阵V_k∈R^N×k包含正交列，Λ_k∈R^k×k是包含矩阵H^TH的k(k≤N)个最大特征值的对角矩阵，T为转置操作，如果正向矩阵H的秩为k，则上式左右完全相等，截断参数k即是用来控制重建精度与运算成本之间的权衡。

于是，对条件协方差矩阵Γ_x＝(κ^-2H^TH+ρI_N)^-1进行近似，并利用Woodbury等式和酉矩阵V_k具有标准正交列的事实，近似的协方差矩阵

可表示为

其中，特征值λ_q(q＝1,...,k)为对角矩阵Λ_k的第q个对角线元素，为了得到均值m_x的近似值，用

替代Γ_x可得

根据以上推论，用于采样的近似目标分布为

利用

的因式分解可对

中进行采样，

因为

是因式分解后对角矩阵，并满足特征值数量k＜＜N，所以可以提供一种计算成本低的方法，来加速高维目标分布

的采样过程。

综上，每次采样关于变量x的闭合解

为

设定采样次数N_samps和老化的样本数N_b，本发明中，取N_samps＝200，N_b为N_samps的10％，初始化参数κ²和参数ν，运用Gibbs采样算法动态构造马尔可夫链，具体抽样步骤如下：

1)通过式(23)计算关于变量x的闭合解

2)利用P-MYULA算法计算第t次迭代中分裂变量z的采样值z^t；

3)计算第t次迭代中变量κ²的采样值

4)计算第t次迭代中变量ν的采样值

5)令t＝t+1并返回步骤1)，直到达到采样次数N_samps；

6)根据

从采样样本中计算出参数x的估计值，即重建的客体密度分布。

以上的实施方式仅是用来说明本发明，对于本技术领域的技术人员来说，在不脱离本发明的基本原理的前提下还可以做出若干改进和润饰，这些改进和润饰都应当视为本发明的保护范围。

Claims (7)

1.一种X射线图像线性重建方法，其特征在于，包括以下步骤：

步骤1：获取X射线投影图像y，引入基于全变差先验的正则项，构造重建的目标函数；

步骤2：引入超先验参数，构造分层贝叶斯模型，建立最大后验概率模型；所述超先验参数为在参数选择之前需要设置的参数，即将参数先验概率分布中的参数在贝叶斯推理过程中看作随机变量进行处理，所述参数先验概率分布中的参数包括均值、方差；所述分层贝叶斯模型为先验概率密度函数含有超先验参数，将所述超先验参数视作随机变量，并服从某一先验分布；

步骤3：运用变量分裂方法，引入分裂变量z，分离数据保真项和全变差TV正则项，得到分裂形式下的联合概率密度分布；所述变量分裂方法是引入分裂变量z，分离数据保真项和全变差TV正则项，然后求解相应的最小化问题；

步骤4：基于Jefferys先验定义超先验变量，得到各变量的条件分布；迭代更新超先验参数，对含全变差TV正则项的分裂变量z的条件分布进行求解；所述Jefferys先验，即非信息先验分布，先验的形式不随着原分布参数形式的变化而变化；

步骤5：计算低秩近似的目标分布得到关于变量x的闭合解，计算采样样本的均值，对变量x进行估计。

2.根据权利要求1所述的一种X射线图像线性重建方法，其特征在于，

在步骤1中，假设投影图像y被加性噪声破坏，则正向问题的线性随机模型为

y＝Hx+e (1)

其中，x∈R^N为需要重建的客体密度分布；y∈R^M为观测的X射线投影图像；H∈R^M×N为数值离散后得到的系统正向矩阵；e∈R^M表示加性噪声或者测量误差，e是一个均值为0、方差为1/μ的高斯随机变量，μ表示噪声参数，R表示实数空间，R^N和R^M×N分别表示N维和M×N维的实数空间；

TV正则项：

其中，||·||₂表示向量的二范数，即向量元素绝对值的平方和再开方，

表示变量x的二维离散梯度算子，在此设定下，关于变量x的后验概率密度分布p(x|y)为

其中，exp[·]表示以e为底数的指数函数，∝表示成比例，基于线性模型式(1)，从观测投影图像y∈R^M中重建出变量x∈R^N是一个病态逆问题，采用TV正则化方法建模，并最小化如下的重建目标函数进行求解：

为重建的客体密度分布；β为正则化参数。

3.根据权利要求1所述的一种X射线图像线性重建方法，其特征在于，

在步骤2中，基于线性模型式(1)，在已知变量x,μ的条件下求投影图像y表示为y|x,μ～N(Hx,μ^-1I)，将投影图像作为模型的一个变量，其中，～表示服从某一分布，N(·)表示高斯分布，I表示单位矩阵，于是可得到似然函数p(y|x,μ)为：

假设变量x的先验分布服从高斯分布，

联合后验概率密度函数为：

4.根据权利要求1所述的一种X射线图像线性重建方法，其特征在于，

在步骤3中，所述变量分裂方法为：

其中，f()代表数据保真项；r()代表一些正则化函数，arg表示对函数求参数，subjectto表示约束条件；

引入正常数ρ以平衡变量x和分裂变量z之间的相似性，然后定义联合概率分布p(x,z|ρ)为

其中，

是一个散度函数，以保证联合概率分布p(x,z|ρ)是合理的，散度函数

是二次的，并满足

根据式(9)可知，关于变量x和分裂变量z的条件概率分布为

结合式(7)，更新后的联合概率密度分布为

5.根据权利要求4所述的一种X射线图像线性重建方法，其特征在于，

在步骤4中,利用方差变量进行模型参数化，即κ²＝μ^-1，τ²＝β^-1，则关于(κ²,τ²)的Jeffreys先验形式为

κ²为尺度不变先验，τ²为数据水平方差尺度先验，定义精度变量ν＝τ²/κ²，进行参数变换后，步骤3中分裂形式下的联合后验概率密度函数(12)更新为

则推导出各个变量的全条件概率密度函数为

其中，Γ(·)表示Gamma分布，有

ν|x,z,κ²,y～1/Γ((N+2)/2,||z||_TV/(2κ²))；

分裂变量z的求解方法为：利用P-MYULA算法对关于分裂变量z的全条件概率分布(18)进行高维采样。

6.根据权利要求5所述的一种X射线图像线性重建方法，其特征在于，在步骤5中，由变量x的全条件密度分布(17)可得，

根据条件分布(19)进行采样得到的样本用x_s＝m_x+Gε表示，其中随机变量ε～N(0,I)，ε为一个均值为0、方差为I高斯随机变量，条件协方差矩阵Γ_x的因式分解G满足Γ_x＝GG^T，通过一个截断的奇异值分解来近似H^TH，即

其中，酉矩阵V_k∈R^N×k包含正交列，Λ_k∈R^k×k是包含矩阵H^TH的k个最大特征值的对角矩阵，T为转置操作，如果正向矩阵H的秩为k，则上式左右完全相等；

对条件协方差矩阵Γ_x＝(κ^-2H^TH+ρI_N)^-1进行近似，并利用Woodbury等式和酉矩阵V_k具有标准正交列的事实，近似的协方差矩阵

表示为

其中，特征值λ_q(q＝1,...,k)为对角矩阵Λ_k的第q个对角线元素，为了得到均值m_x的近似值，用

替代Γ_x可得

同理，采样的近似目标分布为

利用

的因式分解对

中进行采样，

是因式分解后对角矩阵，并满足特征值数量k＜＜N，每次采样关于变量x的闭合解

为

设定采样次数N_samps和老化的样本数N_b，初始化参数κ²和参数ν，运用Gibbs采样算法动态构造马尔可夫链。

7.根据权利要求5所述的一种X射线图像线性重建方法，其特征在于，具体步骤如下：

1)通过式(23)计算关于变量x的闭合解

2)利用P-MYULA算法计算第t次迭代中分裂变量z的采样值z^t；

3)计算第t次迭代中变量κ²的采样值

4)计算第t次迭代中变量ν的采样值

5)令t＝t+1并返回步骤1)，直到达到采样次数N_samps；

6)根据

从采样样本中计算出参数x的估计值，即重建的客体密度分布。