CN110934586B - 一种灰度值矩阵快速分解重建的正则化方法 - Google Patents

Abstract

本发明公开了一种灰度值矩阵快速分解重建的正则化方法，包含如下步骤：(1)根据被测场域获取重建所需的边界测量电压b和灵敏度矩阵A；(2)设置初始化参数；(3)更新辅助变量p^k+1；(4)更新辅助变量q^k+1；(5)更新变量w^k+1；(6)更新变量v^k+1；(7)更新灰度值g^k+1；(8)设置优化参数；(9)更新灰度值u^k+1；(10)更新变量κ^k+1；(11)更新变量m^k+1；(12)判断迭代是否符合迭代终止条件u^k+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(9)步，继续迭代求解；(13)根据最终求解所得灰度值进行成像。本发明能够有效地降低阶梯伪影、提升背景清晰及提高抗噪性。

Description

一种灰度值矩阵快速分解重建的正则化方法

技术领域

本发明属于电阻抗层析成像图像重建方法技术领域，具体涉及一种灰度值矩阵快速分解重建的正则化方法。

背景技术

电学层析成像(Electrical Tomography，ET)是层析成像技术的一种，通过对被测物体施加电流激励获得测量电压数据值，利用特定数学手段逆推被测物体内部的电特性参数分布(电导率/介电系数/复导纳/磁导率)，从而实现场域内部可视化的过程。与磁共振(MRI)，CT等成像技术相比，电学层析成像具有非侵入、无辐射、实时性等优势。例如电阻抗层析成像(EIT)技术，目前已广泛应用于工业过程监测、地球物理地下探测和生物组织诊断等领域。值得注意的是，在医学成像领域上，EIT已被用于乳腺肿瘤的检查、肺通气和脑部成像。

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

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

发明内容

本发明的目的在于提出一种灰度值矩阵快速分解重建的正则化方法，该方法可以有效地降低阶梯伪影、提升背景清晰及提高抗噪性。

本发明为实现上述目的采用如下技术方案，一种灰度值矩阵快速分解重建的正则化方法，其特征在于具体过程为：将电学层析成像看作一个线性不适定问题Ag＝b，其中，A为灵敏度矩阵，b为边界测量电压差值，g为灰度值；设计最小化的目标函数为

式中，

是保真项，λ_TGV是正则化参数，A代表着雅可比矩阵，g代表着灰度值，b为边界测量电压差值，/>

是惩罚项，α₁，α₀是自适应正则化权重因子，

为梯度算子，/>

表示对称梯度算子，v是对称场量；

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

其中，

Ω是图像域集合且/>

是实数域；

由此可推导出

式中σ，τ代表着参数；

为了进一步提高图像空间分辨率，抑制阶梯伪影，针对步骤五迭代的灰度值，本发明提出了一种全新的惩罚项用于图像重建过程。设计了灰度值优化迭代过程：对于优化后灰度值u来说，因为Y^TY＝I，所以可以得到u＝Y^TYu，Y表示灰度值快速分解系数矩阵，Y^T是Y的逆变换，Y^T表示灰度值快速重建系数矩阵。可以利用惩罚项具有沉余性和多分辨率结构的特点，结合L₀范数，对灰度值进行快速分解重建，这样可以稀疏的逼近图像的分段光滑函数，降低阶梯伪影，提高图像空间分辨率。其目标函数定义为：

式中，/>

是保真项，C是可以自适应选择图像不同区域的微分算子，g是灰度值，||λ·κ||₀是惩罚项，λ是正则化参数，κ是优化系数。

为了解决上述L₀范数最小化问题，我们采用拉格朗日算法求解最小化问题。上述式子可转化为

式中m是辅助变量，γ，μ是参数。

通过迭代可得到

由于方程的每一个式子都有一个封闭的形式解，所以L₀极小化问题可以改写为

式中，算子H是一个硬阈值函数。

最终，根据求解优化后灰度值u完成成像。

本发明与现有技术相比具有以下有益效果：本发明的灰度值矩阵快速分解重建的正则化方法在传统TGV正则化的基础上，通过再次迭代灰度值来提高背景清晰度和抑制产生的伪影。另外，为了解决所提出的L₀极小化问题，采用拉格朗日算法进行求解。相比于TV正则化与TGV正则化，本发明提出的灰度值矩阵快速分解重建的正则化方法在降低阶梯伪影和提升图像清晰度等方面有着很好的效果，同时自适应参数的选择增强了参数选择的客观性和简易性。

附图说明

图1为本发明的灰度值矩阵快速分解重建的正则化方法的流程框图。

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

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

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

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

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

具体实施方式

结合附图和实施例对本发明的一种灰度值矩阵快速分解重建的正则化方法加以说明。

本发明提出的一种灰度值矩阵快速分解重建的正则化方法，针对传统正则化算法产生的阶梯伪影问题和背景不清晰问题，以TGV正则化为基础，结合提出的灰度值优化方法，采用L₂范数作为保真度项，L₀范数作为惩罚项，通过自适应选择正则化参数和权重因子来选取最优值，以一种高效、快速的拉格朗日算法来求解本发明提出的目标函数，完成最终的逆问题求解。

如图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正则化为基础，结合本发明提出的灰度值优化方法，通过自适应选择正则化参数和权重因子来选取最优值，以一种高效、快速的拉格朗日算法求解本发明提出的目标函数。

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

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：

其中，A_ij是第j个电极对对第i个电极对的灵敏度系数，/>

分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布。

步骤二：将非线性电阻层析成像的逆问题转化为线性问题；边界测量电压与电导率分布的关系是非线性的f(σ)＝b，式中σ是电导率。对于变化不大的电导分布，f(σ)＝b可将边界测量电压的变化简化为线性形式

式中Δσ是电导率分布的扰动，Δb是电导率的变化引起的边界测量电压的变化。最后导出电阻层析成像逆问题的线性化形式Ag＝b，式中g为所求成像的灰度值。

步骤三：设计目标函数为：

式中，/>

是保真项，λ_TGV是正则化参数，A代表着雅可比矩阵，g代表着灰度值，/>

是惩罚项，α₁，α₀是自适应权重因子，/>

为梯度算子，/>

表示对称梯度算子，v是对称场量。在电阻层析成像图像重建中，通过将目标函数极小化来求解最优灰度值，其模型式表示为：/>

步骤四：推导出模型

的对偶形式为：

式中，

Ω是图像域且/>

是实数域。

步骤五：根据步骤四推导的对偶形式的算法如下

(1)初始化：w＝0,v＝0,

p＝0,/>

q＝0,g₀＝0,τ＝1/L,σ＝1/L

(2)更新辅助变量p^k+1：

(3)更新辅助变量q^k+1：

(4)更新变量w^k+1：

(5)更新变量v^k+1：v^k+1＝v^k+τ(p^k+1+divεq^k+1)。

(6)更新初步灰度值g^k+1：

步骤六：为了进一步提高图像空间分辨率，抑制阶梯伪影，针对步骤五迭代的灰度值，本发明提出了一种全新的惩罚项用于图像重建过程。设计了灰度值优化迭代过程：对于优化后灰度值u来说，因为Y^TY＝I，所以可以得到u＝Y^TYu，Y表示灰度值快速分解系数矩阵，Y^T是Y的逆变换，Y^T表示灰度值快速重建系数矩阵。可以利用惩罚项具有沉余性和多分辨率结构的特点，结合L₀范数，对灰度值进行快速分解重建，这样可以稀疏的逼近图像的分段光滑函数，降低阶梯伪影，提高图像空间分辨率。其目标函数定义为：

式中，/>

是保真项，C是可以自适应选择图像不同区域的微分算子，g是灰度值，||λ·κ||₀是惩罚项，λ是正则化参数，κ是优化系数。

步骤七：利用拉格朗日算法推导出目标函数的模型式

式中，m是辅助变量，μ，γ是设置的参数。

步骤八：根据步骤七推导的模型式的算法如下

(1)初始化：κ＝0，

k＝0，λ＝1e-6，γ＝0.001，μ＝0.001。

(2)更新优化后灰度值u^k+1：u^k+1＝(C^TC+(μ+γ)I)^-1(C^Tg+γu^k+μY^T(κ^k-m^k))。

(3)更新变量κ^k+1：κ^k+1＝H_λ,μ,γ(Yu^k+1+m^k,κ^k)。

(4)更新变量m^k+1：m^k+1＝m^k+Yu^k+1-κ^k+1。

(5)判断迭代是否符合迭代终止条件u^k+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回步骤八第(2)步，继续迭代求解。

步骤九：根据最终求解优化后灰度值u，进行成像。

图3为本发明的实施例选取三种典型的介质模型，TV算法，TGV算法和本发明的方法的图像重建结果示意图。可以看出，三种典型模型中，采用TV正则化方法时，背景中存在明显的阶梯伪影，严重影响了图像恢复的质量。与TV相比，TGV有效地降低了阶梯伪影。然而在边缘保持方面表现不好。相比于前面两种正则化算法，本发明提供的灰度值矩阵快速分解重建的正则化方法很好的抑制了阶梯，并提升了图像背景的清晰程度，同时有着更高的图像空间分辨率。在实际应用中，电学层析成像系统容易受到噪声的影响。为了验证本发明提出的一种灰度值矩阵快速分解重建的正则化方法的抗噪声性能，图4比较了三种方法在不同噪声水平下对模型B的重建图像。可以看出，当噪声水平逐渐增大时，本发明提出的方法在抑制阶梯伪影优于其它两种方法。同时，在电学层析成像中，通常采用图像相对误差(Relative Error，RE)和相关系数(Correlation Coefficient，CC)评价算法来定量图像重建质量，表达式如①、②所示，图像相对误差越小，相关系数越大，表明图像重建质量越好。

式中σ是重建区域的计算电导率，σ^*是实际电导率，t表示像素数，

和/>

表示σ和σ^*的平均值，σ_i和σ_i ^*表示的是σ和σ^*的第i个三角形单元。

图5给出了这三种方法不同噪声水平下模型B的图像相对误差和相关系数，可以看出，本发明提出的一种灰度值矩阵快速分解重建的正则化方法与TV算法，TGV算法相比，具有最低的相对误差和最高的相关系数，能够准确的重建出被测场域内分布，明显地提高了电学层析成像逆问题求解精度和抗噪声能力。

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

Claims (1)

1.一种灰度值矩阵快速分解重建的正则化方法，其特征在于具体过程为：

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

灵敏度矩阵是根据空场状态下的边界测量电压数据得到的，结合灵敏度理论通过公式计算得到灵敏度矩阵，计算公式为：

其中，A_ij是第j个电极对对第i个电极对的灵敏度系数，/>

分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

步骤二：电阻抗层析成像的逆问题是一个高度的非线性问题，需要转化为线性问题进行求解：边界测量电压与电导率分布的关系是非线性的f(σ)＝b，式中σ是电导率，对于变化不大的电导分布，f(σ)＝b可将边界测量电压的变化简化为线性形式

式中Δσ是电导率分布的扰动，Δb是电导率的变化引起的边界测量电压的变化，最后导出电阻层析成像逆问题的线性化形式Ag＝b，式中g为所求成像的灰度值；

步骤三：设计目标函数为：

式中，/>

是保真项，λ_TGV是正则化参数，A代表着雅可比矩阵，g代表着灰度值，b是测量电压，

是惩罚项，α₁，α₀是自适应权重因子，▽为梯度算子，ε(v)＝(▽v+▽v^T)/2表示对称梯度算子，v是对称场量，在电阻抗层析成像图像重建中，通过将目标函数极小化来求解最优灰度值，其模型式表示为：/>

步骤四：推导出模型

的对偶形式为：

式中，

Ω是图像域，/>

是实数域；

步骤五：根据步骤四推导的对偶形式的迭代算法如下：

(1)初始化：w＝0,v＝0,

p＝0,/>

q＝0,g₀＝0,τ＝1/L,σ＝1/L；

(2)更新辅助变量p^k+1：

(3)更新辅助变量q^k+1：

(4)更新变量w^k+1：

(5)更新变量v^k+1：v^k+1＝v^k+τ(p^k+1+divεq^k+1)；

(6)更新灰度值g^k+1：g^k+1＝g^k+τ(div▽p^k+1-A^Tw^k+1)；

步骤六：为了进一步提高图像空间分辨率，抑制阶梯伪影，针对步骤五迭代的灰度值，设计了灰度值优化迭代过程：对于优化后灰度值u来说，因为Y^TY＝I，所以得到u＝Y^TYu，Y表示灰度值快速分解系数矩阵，Y^T是Y的逆变换，Y^T表示灰度值快速重建系数矩阵；利用惩罚项具有沉余性和多分辨率结构的特点，结合L₀范数，对灰度值进行快速分解重建，这样能够稀疏的逼近图像的分段光滑函数，降低阶梯伪影，提高图像空间分辨率，其目标函数定义为：

式中，/>

是保真项，C是可以自适应选择图像不同区域的微分算子，g是灰度值，||λ·κ||₀是惩罚项，λ是正则化参数，κ是优化系数；

步骤七：推导出目标函数的拉格朗日形式：

式中，m是辅助变量，μ，γ是设置的参数；

步骤八：根据步骤七拉格朗日形式的迭代算法如下：

(1)初始化：κ＝0，

k＝0，λ＝1e-6，γ＝0.001，μ＝0.001；

(2)更新优化后灰度值u^k+1：u^k+1＝(C^TC+(μ+γ)I)^-1(C^Tg+γu^k+μ^T(κY^k-m^k))；

(3)更新变量κ^k+11：κ^k+1＝H_λ,μ,γ(Yu^k+1+m^k,κ^k)；

(4)更新变量m^k+1：m^k+1＝m^k+Yu^k+1-κ^k+1；

(5)判断迭代是否符合迭代终止条件u^k+1是否收敛，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回步骤八的第(2)步，继续迭代求解；

步骤九：根据最终求解优化后灰度值u进行成像。

Images