CN110992385B - 抑制伪影与保护边缘的颅内图像重建方法 - Google Patents

Abstract

本发明公开了一种抑制伪影与保护边缘的颅内图像重建方法，包含如下步骤：根据被测场域获取重建所需的相对边界测量值向量b和灵敏度矩阵A；设置初始化参数；更新灰度值g^l+1；更新张量w^l+1；更新辅助变量I P^l+1；更新辅助变量Q^l+1；更新变量更新变量判断迭代是否符合迭代终止条件g^l+1,w^l+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置l＝l+1并跳回更新灰度值步骤，继续迭代求解；根据最终求解所得灰度值进行成像。本发明在降低阶梯伪影和保护边缘等方面有着很好的效果，同时自适应参数的选择增强了参数选择的客观性和简易性。

Description

抑制伪影与保护边缘的颅内图像重建方法

技术领域

本发明属于电学层析成像技术领域，具体涉及利用电学层析成像抑制伪影与保护边缘的颅内图像重建方法以实现图像重建。

背景技术

电学层析成像(Electrical Tomography，ET)是层析成像技术的一种，其通过对被测物体施加电激励，并检测其边界值的变化，利用特定数学手段逆推被测物体内部的电特性参数分布(电导率/介电系数/复导纳/磁导率)，从而得到物体内部的分布情况。与其它层析成像技术相比，电学层析成像具有无辐射、响应速度快、价格低廉等优势。EIT技术已广泛应用于工业过程监测，地球物理地下探测和生物组织诊断。特别是，它在医学成像领域引起了相当大的关注。因此，EIT已被用于乳腺肿瘤的检查、肺通气和脑功能。

电学层析成像逆问题(即图像重建过程)求解具有非线性，通过线性化处理可以将非线性问题转化为线性逆问题求解。由于获取的场域边界电压数量远小于求解场域的像素值，造成逆问题求解的不适定性，所以一般选取正则化方法处理逆问题，例如Tikhonov正则化方法。这种方法能够很好地恢复边缘光滑的物体，而当恢复边缘锐利的物体时则会出现过度的平滑。为了缓解Tikhonov正则化方法带来的过度光滑性，研究了用全变分(Totalvariation regularization method，TV)正则化方法重建尖锐不连续物体的方法。然而，无论原始目标是否连续，都会产生分段不变的图像，称为“阶梯伪影”。为了抑制伪影，对TV正则化方法的改进进行了大量的工作，提出了总广义变分(Total generalized variation，TGV)，在一定程度上缓解了阶梯伪影。

在现有的研究中，这些正则化算法都局限于对噪声敏感的L₂范数保真度项和在一定程度上产生阶梯伪影的惩罚项。这使得成像结果容易受到噪声的影响，并且无法很好的去除阶梯伪影。因此，上述几种正则化方式仍旧无法得到很好的成像效果，而且算法的结构更加复杂，使得很难在实际运用中推广使用。

发明内容

本发明的目的在于提供了一种抑制伪影与保护边缘的颅内图像重建方法，该方法可以有效地降低阶梯伪影、防止边缘过度光滑及提高抗噪性，同时简化了算法结构，提高了算法的求解速度。

本发明为实现上述目的采用如下技术方案，抑制伪影与保护边缘的颅内图像重建方法，其特征在于具体过程为：将电学层析成像看作一个线性不适定问题Ag＝b，其中，A为灵敏度矩阵，b为相对边界测量值向量，g为所求成像灰度值；设计最小化的目标函数为：

式中，是保真项，λ是正则化参数，A代表着雅可比矩阵，g代表着灰度值，是惩罚项，α₁，α₀是自适应正则化权重因子，/>为梯度算子，/>是一个非凸函数，/>Ω是图像域且/> 是微分张量；

在迭代重加权l₁(IRL1)算法基础上将目标函数转化为凸目标函数：

式中，是L₁范数，(g_k,w_k)是上式J_k-1的能量最小值，k是J_k函数的内迭代次数，g_k与w_k是变量g与w内迭代k次的结果；

J_k函数不严格等于原始非凸目标函数，考虑到极小化问题不受常数的影响，在极小化过程中能够忽略，J_k函数的斜率应与原始非凸目标函数在点(g^k,w^k)处的斜率相同；

为了简化上述公式，令上式简化为：

为了解决上述公式的L₁最小化问题，在Chambolle-Pock原始对偶算法的基础上，在最小化目标函数中引入了辅助变量I：p和辅助变量II：q，将最小化目标函数的对偶形式根据Legendre-Fenchel对偶形式导出：

其中，

由此推导出

由于求解原始对偶算法时收敛速度，故在迭代过程中加入了松弛算法，以便于快速迭代出所需灰度值，提高整个算法的效率：

式中，θ代表松弛参数；

重建算法的具体过程为：(1)根据被测场域获取重建所需的相对边界测量值向量b和灵敏度矩阵A；边界测量值的获取是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b为不含内含物的空场边界测量电压向量b_meas1和含有内含物的有物场的边界测量电压向量b_meas2之差；(2)设置初始化参数：调节最优权重因子的正则化参数λ，调节最优非凸函数指数t，最大迭代次数l_max，最大内迭代次数k_max，自适应正则化权重因子α₁，α₀；(3)更新灰度值g^l+1；(4)更新张量ω^l+1；(5)更新辅助变量I：P^l+1；(6)更新辅助变量II：Q^l+1；(7)更新变量(8)更新变量/>(9)判断迭代是否符合迭代终止条件g^l+1,w^l+1是否收敛，若是，则迭代终止，进行下一步操作；若否，设置l＝l+1并跳回第(3)步，继续迭代求解；(10)根据最终求解所得灰度值进行成像。

进一步优选，所述重建算法的具体步骤为：

(1)针对三种典型模型分别获取各自重建所需的边界测量值b和灵敏度矩阵A：边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值，104个独立测量值；b为不含内含物的空场边界电压b_meas1和含有内含物的有物场的边界测量电压b_meas2之差即右端项相对边界测量值b＝b_meas2-b_meas1；

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理计算灵敏度矩阵，计算公式为：其中A_ij是第j个电极对对第i个电极对的灵敏度系数，/>分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：λ是正则化参数，调节最优非凸函数指数t＝0.7，初始迭代次数l＝0，最大迭代次数l_max＝1000，初始内迭代次数k＝0，最大内迭代次数k_max＝200，自适应正则化权重因子α₁＝2.5和α₀＝0.3，初始灰度值g⁰＝A^Tb；

(3)更新灰度值g^l+1：P^l＝m_kp^l；

(4)更新张量w^l+1：Q^l＝n_kq^l；

(5)更新辅助变量IP^l+1：

(6)更新辅助变量IIQ^l+1：

(7)更新变量

(8)更新变量

(9)判断迭代是否符合迭代终止条件g^l+1,w^l+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置l＝l+1并跳回第(3)步，继续迭代求解；

(10)根据最终求解所得灰度值进行成像。

本发明与现有技术相比具有以下有益效果：本发明提出的抑制伪影与保护边缘的颅内图像重建方法在传统TGV正则化的基础上，通过在惩罚项中加入非凸函数来抑制伪影和保护边缘。另外，为了解决所提出的非凸极小化问题，采用迭代重加权L₁(IRL1)算法将非凸目标函数转换为凸目标函数，并用Chambolle-Pock原始对偶算法进行求解。相比于TV正则化与TGV正则化，本发明提出的抑制伪影与保护边缘的颅内图像重建方法在降低阶梯伪影和保护边缘等方面有着很好的效果，同时自适应参数的选择增强了参数选择的客观性和简易性。

附图说明

图1为本发明的抑制伪影与保护边缘的颅内图像重建方法的流程框图。

图2为本发明的电阻层析成像系统环形单截面被测场域，激励电流和测量电压的模式以及电极分布。

图3为本发明的实施例在选取三种模型的真实分布时，TV正则化算法，TGV正则化算法和抑制伪影与保护边缘的颅内图像重建方法的图像重建结果图。

图4为三种方法在不同噪声水平下对模型2的重建图像。

图5为三种方法在不同噪声水平下模型2的图像相对误差和相关系数。

图中：1-被测场域，2-电极，3-激励电流，4-测量电压。

具体实施方式

结合附图和实施例对本发明的一种抑制伪影与保护边缘的颅内图像重建方法加以说明。

本发明提出的一种抑制伪影与保护边缘的颅内图像重建方法，针对传统正则化算法产生的阶梯伪影和边缘过度光滑问题，以TGV惩罚项和非凸正则化为基础，结合采用L₂范数作为保真度项，L₁范数作为惩罚项，通过自适应选择正则化参数和权重因子来选取最优值，以一种高效、快速的迭代重加权L₁(IRL1)算法和Chambolle-Pock原始对偶算法来求解本发明提出的目标函数，完成最终的逆问题求解。

如图1所示，为本发明的抑制伪影与保护边缘的颅内图像重建方法流程图。

如图2所示，为电学层析成像之一的生物电阻抗层析成像系统环形单截面被测场域1，激励电流3和测量电压4的模式以及电极2分布，采用16电极2均匀分布在被测场域1外壁。

选取三种典型的介质模型为实施例，被测场域内物体真实分布如图3左侧一竖列所示，其他三列从左到右分别表示为TV正则化算法，TGV正则化算法和抑制伪影与保护边缘的颅内图像重建方法。为了较好地体现本发明中一种抑制伪影与保护边缘的颅内图像重建方法与其它两种方法的不同，分别给出了三种典型模型在这三种正则化算法下的求解结果。

具体实施过程包括如下步骤：抑制伪影与保护边缘的颅内图像重建方法，该方法将电学层析成像看作一个线性不适定问题Ag＝b。其中，A为灵敏度矩阵，b为相对边界测量值向量，g代表着灰度值。逆问题可以用最小二乘优化形式的目标函数表示：其中，f(g)为目标函数。正则化方法为提高解的稳定性提供了一种新的选择，其形式可以描述为：/>其中，λ是正则化参数。

TV正则化算法在保持锐利边缘方面表现出良好的效果。其一般形式可以表示为：其中，/>表示TV正则化算法中的惩罚项。在边缘保持方面，TV正则化方法有着明显的优势。然而TV正则化分段逼近，在图像重建中可以观察到明显的阶梯伪影。

为了克服传统TV正则化方法的不足，提高重建图像的质量，TGV正则化被提了出来，其一般形式可以表示为：其中/>表示TGV正则化算法中的惩罚项。在一定程度上抑制了阶梯伪影，但是在边缘保持方面表现不好。

本发明针对传统正则化算法产生的阶梯伪影和边缘过度光滑问题，以TGV惩罚项和非凸正则化为基础，通过自适应选择正则化参数和权重因子来选取最优值，以一种高效、快速的迭代重加权L₁(IRL1)算法和Chambolle-Pock原始对偶算法来求解本发明提出的目标函数。其最小化的目标函数为：

式中A代表着雅可比矩阵，g代表着灰度值，λ是自适应正则化参数，是惩罚项，α₁，α₀是自适应正则化权重因子。/>是为梯度算子，/>是一个非凸函数，/>Ω是图像域且/> 是一阶微分张量。

在迭代重加权L₁(IRL1)算法基础上将目标函数转化为凸目标函数

J_k函数不严格等于原始非凸目标函数，但考虑到极小化问题不受常数的影响，在极小化过程中是可以忽略的。然而，J_k函数的斜率应与原始非凸目标函数在点(g^k,w^k)处的斜率相同。

为了简化上述式子，令上式可以简化为：

为了解决上述式子的L₁最小化问题，在Chambolle-Pock原始对偶算法的基础上，在最小化目标函数中引入了辅助变量I：p和辅助变量II：q，将最小化目标函数的对偶形式根据Legendre-Fenchel对偶形式导出并求解：

式中

由此可以推导出

重建算法的具体步骤为：

(1)针对三种典型模型分别获取各自重建所需的边界测量值b和灵敏度矩阵A：边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极(如图2所示)，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值，104个独立测量值；b为不含内含物的空场边界电压b_meas1和含有内含物的有物场的边界测量电压b_meas2之差即右端项相对边界测量值b＝b_meas2-b_meas1；

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：其中A_ij是第j个电极对对第i个电极对的灵敏度系数，/>分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：λ是正则化参数，调节最优非凸函数指数t＝0.7，初始迭代次数l＝0，最大迭代次数l_max＝1000，初始内迭代次数k＝0，最大内迭代次数k_max＝200，自适应正则化权重因子α₁＝2.5和α₀＝0.3，初始灰度值g⁰＝A^Tb；

(3)更新灰度值g^l+1：P^l＝m_kp^l；

(4)更新张量w^l+1：Q^l＝n_kq^l；

(5)更新辅助变量IP^l+1：

(6)更新辅助变量IIQ^l+1：

(7)更新变量

(8)更新变量

(9)判断迭代是否符合迭代终止条件g^l+1,w^l+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置l＝l+1并跳回第(3)步，继续迭代求解；

(10)根据最终求解所得灰度值进行成像。

图3为本发明的实施例，选取三种典型的介质模型，TV算法，TGV算法和抑制伪影和保护边缘的颅内图像重建方法的图像重建结果图。可以看出，三种典型模型中，采用TV正则化方法时，背景中存在明显的阶梯伪影，严重影响了图像恢复的质量。与TV相比，TGV有效地降低了阶梯伪影。然而在边缘保持方面表现不好。相比于前面两种正则化算法，一种抑制伪影和保护边缘的颅内图像重建方法的重建图像的阶梯伪影得到了很好的抑制，并防止重建图像边缘过度光滑，并且有着更高的图像分辨率，明显的提高了图像重建的质量。在实际应用中，电学层析成像系统容易受到噪声的影响。

为了验证本发明提出的电学层析成像伪影抑制图像重建方法的抗噪声性能，图4比较了三种方法在不同噪声水平下对模型B的重建图像。可以看出，当噪声水平逐渐增大时，本发明提出的方法在抑制阶梯伪影和保持锐利边缘方面优于其它两种方法。同时，在电学层析成像中，通常采用图像相对误差(Relative Error，RE)和相关系数(CorrelationCoefficient，CC)评价算法来定量评价图像重建质量，表达式如①、②所示，图像相对误差越小，相关系数越大，表明图像重建质量越好。

式中σ是重建区域的计算电导率，σ^*是实际电导率，t表示像素数，和/>表示σ和σ^*的平均值，σ_i和σ_i ^*表示的是σ和σ^*的第i个三角形单元。

图5给出了这三种方法不同噪声水平下模型B的图像相对误差和相关系数，可以看出，本发明提出的一种抑制伪影和保护边缘的颅内图像重建方法与TV算法，TGV算法相比，具有最低的相对误差和最高的相关系数，能够准确的描绘出被测场域内部分布，明显地提高了电学层析成像逆问题求解精度和抗噪声能力。

以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.抑制伪影与保护边缘的颅内图像重建方法，其特征在于具体过程为：将电学层析成像看作一个线性不适定问题Ag＝b，其中，A为灵敏度矩阵，b为相对边界测量值向量，g为所求成像灰度值；设计最小化的目标函数为：

式中，是保真项，λ是正则化参数，A代表着雅可比矩阵，g代表着灰度值，是惩罚项，α₁，α₀是自适应正则化权重因子，/>为梯度算子，/>是一个非凸函数，/>Ω是图像域且/> 是微分张量；

在迭代重加权l₁(IRL1)算法基础上将目标函数转化为凸目标函数：

式中，是L₁范数，(g_k,w_k)是上式J_k-1的能量最小值，k是J_k函数的内迭代次数，g_k与w_k是变量g与w内迭代k次的结果；

J_k函数不严格等于原始非凸目标函数，考虑到极小化问题不受常数的影响，在极小化过程中能够忽略，J_k函数的斜率应与原始非凸目标函数在点(g^k,w^k)处的斜率相同；

为了简化上述公式，令上式简化为：

为了解决上述公式的L₁最小化问题，在Chambolle-Pock原始对偶算法的基础上，在最小化目标函数中引入了辅助变量I：p和辅助变量II：q，将最小化目标函数的对偶形式根据Legendre-Fenchel对偶形式导出：

其中，

由此推导出

由于求解原始对偶算法时收敛速度，故在迭代过程中加入了松弛算法，以便于快速迭代出所需灰度值，提高整个算法的效率：

式中，θ代表松弛参数；

重建算法的具体过程为：(1)根据被测场域获取重建所需的相对边界测量值向量b和灵敏度矩阵A；边界测量值的获取是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b为不含内含物的空场边界测量电压向量b_meas1和含有内含物的有物场的边界测量电压向量b_meas2之差；(2)设置初始化参数：调节最优权重因子的正则化参数λ，调节最优非凸函数指数t，最大迭代次数l_max，最大内迭代次数k_max，自适应正则化权重因子α₁，α₀；(3)更新灰度值g^l+1；(4)更新张量ω^l+1；(5)更新辅助变量I：P^l+1；(6)更新辅助变量II：Q^l+1；(7)更新变量(8)更新变量/>(9)判断迭代是否符合迭代终止条件g^l+1,w^l+1是否收敛，若是，则迭代终止，进行下一步操作；若否，设置l＝l+1并跳回第(3)步，继续迭代求解；(10)根据最终求解所得灰度值进行成像。

2.根据权利要求1所述的抑制伪影与保护边缘的颅内图像重建方法，其特征在于所述重建算法的具体步骤为：

(1)针对三种典型模型分别获取各自重建所需的边界测量值b和灵敏度矩阵A：边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值，104个独立测量值；b为不含内含物的空场边界电压b_meas1和含有内含物的有物场的边界测量电压b_meas2之差即右端项相对边界测量值b＝b_meas2-b_meas1；

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理计算灵敏度矩阵，计算公式为：其中A_ij是第j个电极对对第i个电极对的灵敏度系数，分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：λ是正则化参数，调节最优非凸函数指数t＝0.7，初始迭代次数l＝0，最大迭代次数l_max＝1000，初始内迭代次数k＝0，最大内迭代次数k_max＝200，自适应正则化权重因子α₁＝2.5和α₀＝0.3，初始灰度值g⁰＝A^Tb；

(3)更新灰度值g^l+1：P^l＝m_kp^l；

(4)更新张量w^l+1：Q^l＝n_kq^l；

(5)更新辅助变量IP^l+1：

(6)更新辅助变量IIQ^l+1：

(7)更新变量

(8)更新变量

(9)判断迭代是否符合迭代终止条件g^l+1,w^l+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置l＝l+1并跳回第(3)步，继续迭代求解；

(10)根据最终求解所得灰度值进行成像。