CN111062999B - 有效保留锐利边缘的生物医学电阻抗层析成像方法 - Google Patents

Abstract

本发明公开了有效保留锐利边缘的生物医学电阻抗层析成像方法，包含如下步骤：(1)根据被测场域获取重建所需的相对边界测量电压b和雅克比矩阵A；(2)将非线性的问题转化为线性问题；(3)设置初始化参数；(4)更新权重因子ω；(5)更新非凸参数π₁，π₂；(6)更新正则化参数λ；(7)更新辅助变量v,w,z；(8)更新电导率分布g；(9)判断迭代是否符合迭代终止条件或者k＞k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1继续迭代求解；(10)根据最终求解所得电导率值进行成像。本发明有效解决混合全变分正则化方法不能有效保留锐利边缘的问题，采用自适应方法的选取正则化参数和权重因子，进一步提高成像的质量。

Description

有效保留锐利边缘的生物医学电阻抗层析成像方法

技术领域

本发明属于电学层析成像技术领域，具体涉及利用有效保留锐利边缘的生物医学电阻抗层析成像方法以实现图像的重建。

背景技术

电学层析成像技术(Electrical Tomography，ET)出现于20世纪80年代后期，是一种基于被测物体电学特性进行图像重建的过程层析成像技术，它通过边界测量值得出被测区域内介质的分布信息，进而对电学特性的分布信息进行成像。电学层析成像被广泛应用在工业过程检测和生物医学成像领域。电学层析成像技术主要包括电阻层析成像(Electrical Resistance Tomography，ERT)、电阻抗层析成像(Electrical ImpedanceTomography，EIT)、电磁层析成像(Electrical Magnetic Tomography，EMT)和电容层析成像(Electrical Capacitance Tomography，ECT)。其中，电阻抗层析成像在生物医学成像领域的临床应用上拥有巨大的潜力，具有非侵入、无辐射、实时性、低成本和便携带等优点。电阻抗层析成像的图像重建是一个严重不适定的逆向问题，而且生物体的结构复杂和外界影响因素较多。因此，电阻抗层析成像在生物医学成像上的应用对重建方法有较高的要求，需要一种稳定、高效、高分辨率的方法进行逆向问题的求解。正则化方法是目前解决不适定问题比较常用方法，通过在目标函数数据保真项的基础上引入惩罚项的方法来稳定解。全变分正则化方法可以有效保留边界的不连续信息，允许重建锐利的边缘；例如G Gonzale等人2018年发表于《科学与工程中的逆向问题》(Inverse Problems in Science andEngineering)第74卷，第564-576页，题为电阻抗层析成像中的各向同性和各向异性全变分正则化(Isotropic and anisotropic total variation regularization in electricalimpedance tomography。然而，全变分正则化方法在重建图像的平滑区域会出现严重的阶梯效应，降低了成像的分辨率。为了解决阶梯效应的问题，对该方法进行了大量改进的工作，其中混合全变分正则化方法在克服阶梯效应和保留边缘信息方面提供了一种较为有效的折中方案；例如T Adam等人2019年发表于《多维系统与信号处理》(MultidimensionalSystems and Signal Processing)第30卷，第503-527页，题为基于重叠群稀疏的高阶非凸总变分的图像去噪(Image denoising using combined higher order non-convex totalvariation with overlapping group sparsity)。混合全变分正则化方法在一阶梯度惩罚项的基础上引入二阶梯度惩罚项，并通过权重因子控制两项之间的权重；希望在抑制阶梯效应与保留边缘信息之间建立折中，但该方法不可避免的牺牲了全变分正则化方法在保留边缘信息方面的性能。为了弥补混合全变分正则化方法保边性能不理想的缺陷，本发明提出了有效保留锐利边缘的生物医学电阻抗层析成像方法，该方法利用非凸函数有利于保留边缘信息的特性，在混合全变分正则化方法惩罚项的基础上引入非凸函数，进一步增强在保留边缘信息方面的性能。

发明内容

本发明的目的在于提出了一种有效保留锐利边缘的生物医学电阻抗层析成像方法，用于有效解决混合全变分正则化方法不能有效保留锐利边缘的问题；该方法在混合全变分正则化方法惩罚项的基础上引入非凸函数，不仅能够有效抑制阶梯效应，还能够很好地保留重建图像的边缘信息；采用自适应方法的选取正则化参数和权重因子，进一步提高成像的质量。

本发明为实现上述目的采用如下技术方案，有效保留锐利边缘的生物医学电阻抗层析成像方法，其特征在于具体过程为：将电阻抗层析成像的非线性逆向问题转化为线性问题Ag＝b，式中，A为灵敏度矩阵，g为真实电导率分布，b为相对边界测量电压；所设计目标函数的优化模型为：

式中，为优化得到的电导率分布，/>为数据保真项又称为最小二乘项，λ为正则化参数用于权衡数据保真项与惩罚项又称为正则项之间的权重；/>为一阶非凸惩罚项，/>为二阶非凸惩罚项，▽为梯度算子，ω为权重因子用于权衡一阶非凸惩罚项与二阶非凸惩罚项之间的权重，0＜π₁＜1和0＜π₂＜1为非凸函数的非凸参数；

非凸目标函数的求解首先需要利用迭代重加权最小二乘法将其变换为凸逼近模型的形式，凸逼近模型的优化模型表示为：

式中，为了避免/>和/>的分母为0，将φ的大小设置为10^-6；

在分裂变量法的基础上引入了3个辅助变量v、w和z，将优化模型的无约束问题转化为约束问题表示为：

s.t.z＝Ag-b,v＝▽g,w＝▽²g

为了用交替方向乘子迭代算法求上式，将上式表示为增广拉格朗日函数：

式中，T表示矩阵的转置算子，α,β和χ是增广拉格朗日乘子，δ₁,δ₂,δ₃>0是平衡各项之间权重的惩罚参数；

重建算法的具体过程为：(1)获取重建所需的相对边界测量电压b和灵敏度矩阵A；采用相邻电流激励和相邻电压测量且激励电极不测量的模式，并在循环激励、循环测量下采集测量电极上的边界电压，相对边界测量电压b为含有内含物的有物场边界测量电压b₁之差与不含内含物的空场边界测量电压b₂；(2)将图像重建的非线性问题转化为线性问题；(3)设置初始化参数；(4)更新权重因子ω；(5)更新非凸参数π₁、π₂；(6)更新正则化参数λ；(7)更新辅助变量v,w,z；(8)更新电导率分布g；(9)更新增广拉格朗日乘子α,β,和χ；(10)判断迭代是否符合迭代终止条件或者k＞k_max，若是，则迭代终止，进行下一步操作；若否，设置k＝k+1继续迭代求解；(11)根据最终求解所得电导率分布/>进行成像。

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

(1)分别获取三种被测模型各自重建所需的相对边界测量电压和雅可比矩阵A；在被测生物体表面均匀分布16个电极，采用相邻模式的电流激励和电压测量，通过采集循环激励和循环测量下的边界电压，共获得208组测量值；相对边界测量电压b为含有内含物的满场边界测量电压b₂与不含内含物的空场边界测量电压b₁之差，即b＝b₂-b₁；

基于有限元仿真的雅可比矩阵的计算公式为：式中，A_ij是第j个电极对对第i个电极对的灵敏度系数，/>分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)将图像重建的非线性问题转化为线性问题：边界测量电压与电导率分布的关系是非线性的，对于变化不大的电导率分布，可将边界测量电压的变化简化为线性形式式中，Δg是电导率分布的扰动，ΔU是电导率的变化引起的边界测量电压的变化，将以上线性化形式表示为b＝Ag；

(3)设置初始化参数：最小迭代阈值ε＝10^-6，初始迭代次数k＝0，最大迭代次数k_max＝200，惩罚参数δ₁＝200，δ₂＝δ₃＝30，增广拉格朗日乘子α⁰＝0，β⁰＝0，χ＝0初始灰度值g⁰＝A^Tb；

(4)更新权重因子ω：根据上一次优化得到的值自动选择合适的权重因子，

(5)更新非凸参数π₁，π₂：根据上次迭代结果值的一阶梯度和二阶梯度自动选择合适的非凸参数，/>

(6)更新正则化参数λ：根据上次迭代的值的最大值和最小值为先验经验自适应的选取合适的正则化参数，/>式中，r代表有限元剖分网格的平均长度，t为调节正则化参数的因子0＜t＜100；

(7)更新辅助变量v：

(8)更新辅助变量w：

(9)更新辅助变量z：

(10)更新电导率分布g：

为了提高计算效率和提高解的稳定性，采用快速傅里叶变换求解子问题：

式中，上标T是转置算子，▽^T▽,(▽²)^T▽²和A^TA是具有循环块结构的块循环矩阵；

(11)更新增广拉格朗日乘子α,β,和χ：

(12)判断迭代是否符合迭代终止条件或者k＞k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1继续迭代求解；

(13)根据最终求解所得电导率值进行成像。

本发明具有以下有益效果：本发明的有效保留锐利边缘的生物医学电阻抗层析成像方法，通过在混合全变分正则化方法的惩罚项中引入非凸函数，进一步提高混合全变分正则化方法保留边缘信息的性能。该方法通过自适应方法选择权重因子来控制一阶非凸惩罚项与二阶非凸惩罚项之间的权重和自适应的选择正则化参数调节数据保真项和惩罚项之间的权重；自适应参数的选择增强了参数选择的客观性，有利于提高成像的分辨率。为了解决所提出的极小化问题，研究了一种高效、快速的交替方向乘子法，该方法简化了算法的复杂度，提高了算法的运行速度，增强了算法的适用性。相比于全变分正则化方法和混合全变分正则化方法，本发明提出的有效保留锐利边缘的生物医学电阻抗层析成像方法在降低阶梯效应、增强保留边缘信息性能和提高成像分辨率等方面有着很好的效果。

附图说明

图1为本发明有效保留锐利边缘的生物医学电阻抗层析成像方法图像重建的流程框图；

图2为选取的两类边缘的模型分别用全变分正则化方法、混合全变分正则化方法和有效保留锐利边缘的生物医学电阻抗层析成像方法图像重建的结果图；

图3为选取的混合模型用全变分正则化方法、混合全变分正则化方法和有效保留锐利边缘的生物医学电阻抗层析成像方法在不同噪声水平下图像重建的结果图；

图4为混合模型重建图像的相对误差(relative error，re)和相关系数(correlation coefficient，cc)；

图5为人体肺部模型分别用全变分正则化方法、混合全变分正则化方法和有效保留锐利边缘的生物医学电阻抗层析成像方法图像重建的结果图。

具体实施方式

结合附图详细对本发明的有效保留锐利边缘的生物医学电阻抗层析成像方法加以说明。

有效保留锐利边缘的生物医学电阻抗层析成像方法，针对利用混合全变分正则化方法不能有效保留重建图像边缘信息的问题。在混合全变分正则化方法的惩罚项的基础上引入非凸函数，进一步提高混合全变分正则化方法保留边缘信息的性能。该方法通过自适应的方法选择权重因子来控制一阶非凸惩罚项与二阶非凸惩罚项之间的权重和自适应的选择正则化参数调节数据保真项和惩罚项之间的权重。为了解决所提出的极小化问题，研究了一种高效、快速的交替方向乘子法对所提出的新方法进行求解，并完成最终图像重建。

如图1所示，为本发明的有效保留锐利边缘的生物医学电阻抗层析成像方法图像重建的流程图。

如图2所示，本发明对两种不同边缘的模型，具有平滑边缘的模型a和具有锐利边缘的模型b进行了研究。图中第一列为重建的真实模型，右边三列依次为全变分正则化方法、混合全变差正则化方法和有效保留锐利边缘的生物医学电阻抗层析成像方法重建的电导率分布的图像。结果表明，采用全变分正则化方法重建的图像会出现较为严重的阶梯效应；而混合全变分正则化方法重建的图像，则有效抑制了阶梯效应，但重建图像的边缘信息未能得到较好的保留，重建目标的轮廓在一定程度上的发生了形变。本发明提出的新方法既能够抑制平滑区域产生的阶梯效应，又能够有效保留锐利边缘的细节信息，与前两种正则化方法相比提高了重建图像的分辨率和整体成像的质量。

电阻抗层析成像的逆向问题是严重不适定性的。目前通常采用正则化方法来解决这类问题，通过在目标函数数据保真项的基础上增加惩罚项，提高解的稳定性。正则化方法的一般形式可以描述为：式中，λ是正则化参数控制着数据保真项和惩罚项R(g)之间权重。

全变分正则化方法是一种经典的正则化方法，在保留边缘信息方面有较好的性能，目标函数的优化模型可以表示为：然而，该方法仅能对分段常数函数进行重建，所以在重建图像的平滑区域会产生明显的阶梯效应。

混合全变分正则方法是在全变分正则化方法的一阶梯度惩罚项的基础上引入二阶梯度惩罚项，两项共同构成目标函数的惩罚项，来抑制在全变分正则化方法重建的图像中产生的阶梯效应，并通过权重因子控制一阶梯度惩罚项和二阶梯度惩罚项之间的权重，目标函数的优化模型可以表示为：该方法虽然有效的抑制了阶梯效应，但在一定上牺牲了全变分正则化方法在保留边缘信息方面的性能。

为了解决电阻抗层析成像在生物医学成像中重建图像分辨率不高的问题，本发明针对混合全变分正则化方法在进行图像重建过程中可以有效抑制阶梯效应，但保边性能较差的问题；提出了有效保留锐利边缘的生物医学电阻抗层析成像方法，该方法利用非凸函数有利于保留边缘信息的特性，在混合全变分正则化方法正则化项的基础上加入非凸函数，并通过自适应的方法对正则化参数、权重因子和非凸参数进行选取，以达到最理想的成像效果。提出一种交替方向乘子法对所提出新方法的目标函数进行求解。设计目标函数的优化模型为：

式中，为优化得到的电导率；/>为数据保真项，g为真实电导率，A为雅可比矩阵，b相对边界测量电压，λ为正则化参数用于权衡数据保真项与惩罚项之间的权重；/>为非凸一阶梯度惩罚项，/>为非凸二阶梯度惩罚项，▽为梯度算子，ω为权重因子，用于权衡一阶非凸惩罚项与二阶非凸惩罚项之间的权重，0＜π₁＜1和0＜π₂＜1为非凸参数。

非凸目标函数利用迭代重加权最小二乘法变换为凸逼近函数。凸逼近模型的优化模型表示为：

其中，为了避免/>和/>的分母为0，将φ的大小设置为10^-6。

在分裂变量法的基础上引入了3个辅助变量v、w和z，将优化模型的无约束问题转化为约束问题表示为：

s.t.z＝Ag-b,v＝▽g,w＝▽²g

为了用交替方向乘子迭代算法求上式，将上式表示为增广拉格朗日函数：

式中，T表示矩阵的转置算子，α,β和χ是增广拉格朗日乘子，δ₁,δ₂,δ₃>0是平衡各项之间权重的惩罚参数。

利用交替方向乘子法求解新方法目标函数，并进行图像重建的具体步骤如下：

(1)分别获取三种被测模型各自重建所需的相对边界测量电压b和雅可比矩阵A；在被测生物体表面均匀分布16个电极，采用相邻模式的电流激励和电压测量，通过采集循环激励和循环测量下的边界电压，共获得208组测量值；相对边界测量电压为含有内含物的满场边界测量电压b₁与不含内含物的空场边界测量电压b₂之差，即：b＝b₁-b₂。

基于有限元仿真的雅可比矩阵的计算公式为：式中，A_ij是第j个电极对对第i个电极对的灵敏度系数，/>分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)将图像重建的非线性问题转化为线性问题:边界测量电压与电导率分布的关系是非线性的。对于变化不大的电导率分布，可将边界测量电压的变化简化为线性形式式中，Δg是电导率分布的扰动，ΔU是电导率的变化引起的边界测量电压的变化。最后为了方便起见，将以上线性化形式表示为b＝Ag；

(3)设置初始化参数：最小迭代阈值ε＝10^-6，初始迭代次数k＝0，最大迭代次数k_max＝200，惩罚参数δ₁＝200，δ₂＝δ₃＝30，增广拉格朗日乘子α⁰＝0,β⁰＝0,χ＝0初始灰度值g⁰＝A^Tb；

(4)更新权重因子ω：根据上一次优化得到的值自动选择合适的权重因子，

(5)更新非凸参数π₁，π₂：根据上次迭代结果值的一阶梯度和二阶梯度自动选择合适的非凸参数，/>

(6)更新正则化参数λ：根据上次迭代值的最大值/>和最小值/>为先验经验自适应的选取合适的正则化参数，/>式中，r代表有限元剖分网格的平均长度，t为调节正则化参数的因子0＜t＜100；

(7)更新辅助变量v：

(8)更新辅助变量w：

(9)更新辅助变量z：

(10)更新电导率分布g：

为了提高计算效率和提高解的稳定性，采用快速傅里叶变换求解子问题，

式中，上标T是转置算子，▽^T▽,(▽²)^T▽²和A^TA是具有循环块结构的块循环矩阵；

(11)更新增广拉格朗日乘子α,β,和χ：

(12)判断迭代是否符合迭代终止条件或者k＞k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1继续迭代求解；

(13)根据最终求解所得电导率值进行成像。

图3比较了三种方法在不同噪声水平下对模型c的重建图像。可以看出，当噪声水平逐渐增大时，本发明提出的有效保留锐利边缘的生物医学电阻抗层析成像方法在抑制阶梯伪影和保持锐利边缘信息方面优于其它两种方法。在电学层析成像中，通常采用图像相对误差(relative error，re)和相关系数(correlation coefficient，cc)评价算法来定量重建图像的质量，图像相对误差越小，相关系数越大，表明图像重建质量越好。相对误差re和相关系数cc的表达式为：

式中，n表示像素数，和/>表示/>和g的平均值，/>和g_i表示的是/>和g的第i个三角形单元。

图4给出了这三种方法不同噪声水平下模型c的图像相对误差和相关系数，可以看出，本发明提出的有效保留锐利边缘的生物医学电阻抗层析成像方法，与全变分正则化方法和混合全变分正则化方法相比，具有最低的相对误差和最高的相关系数，进一步表明所提出的新方法在重建图像时的优越性。

图5为人体肺部模型d分别用全变分正则化方法、混合全变分正则化方法和有效保留锐利边缘的生物医学电阻抗层析成像方法重建的结果图；可以看出在对人体肺部进行重建时，由本发明所提出的新方法重建的图像背景清晰，轮廓分明，更能接近真实的肺部结构。

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

Claims (2)

1.有效保留锐利边缘的生物医学电阻抗层析成像方法，其特征在于具体过程为：将电阻抗层析成像的非线性逆向问题转化为线性问题Ag＝b，式中，A为灵敏度矩阵，g为真实电导率分布，b为相对边界测量电压；所设计目标函数的优化模型为：

式中，为优化得到的电导率分布，/>为数据保真项又称为最小二乘项，λ为正则化参数用于权衡数据保真项与惩罚项又称为正则项之间的权重；/>为一阶非凸惩罚项，为二阶非凸惩罚项，/>为梯度算子，ω为权重因子用于权衡一阶非凸惩罚项与二阶非凸惩罚项之间的权重，0＜π₁＜1和0＜π₂＜1为非凸函数的非凸参数；

非凸目标函数的求解首先需要利用迭代重加权最小二乘法将其变换为凸逼近模型的形式，凸逼近模型的优化模型表示为：

式中，为了避免/>和/>的分母为0，将φ的大小设置为10^-6；

在分裂变量法的基础上引入了3个辅助变量v、w和z，将优化模型的无约束问题转化为约束问题表示为：

s.t.z＝Ag-b,

为了用交替方向乘子迭代算法求上式，将上式表示为增广拉格朗日函数：

式中，T表示矩阵的转置算子，α,β和χ是增广拉格朗日乘子，δ₁,δ₂,δ₃>0是平衡各项之间权重的惩罚参数；

重建算法的具体过程为：(1)获取重建所需的相对边界测量电压b和灵敏度矩阵A；采用相邻电流激励和相邻电压测量且激励电极不测量的模式，并在循环激励、循环测量下采集测量电极上的边界电压，相对边界测量电压b为含有内含物的有物场边界测量电压b₁之差与不含内含物的空场边界测量电压b₂；(2)将图像重建的非线性问题转化为线性问题；(3)设置初始化参数；(4)更新权重因子ω；(5)更新非凸参数π₁、π₂；(6)更新正则化参数λ；(7)更新辅助变量v,w,z；(8)更新电导率分布g；(9)更新增广拉格朗日乘子α,β,和χ；(10)判断迭代是否符合迭代终止条件或者k＞k_max，若是，则迭代终止，进行下一步操作；若否，设置k＝k+1继续迭代求解；(11)根据最终求解所得电导率分布/>进行成像。

2.根据权利要求1所述的有效保留锐利边缘的生物医学电阻抗层析成像方法，其特征在于所述重建算法的具体步骤为：

(1)分别获取三种被测模型各自重建所需的相对边界测量电压和雅可比矩阵A；在被测生物体表面均匀分布16个电极，采用相邻模式的电流激励和电压测量，通过采集循环激励和循环测量下的边界电压，共获得208组测量值；相对边界测量电压b为含有内含物的满场边界测量电压b₂与不含内含物的空场边界测量电压b₁之差，即b＝b₂-b₁；

基于有限元仿真的雅可比矩阵的计算公式为：式中，A_ij是第j个电极对对第i个电极对的灵敏度系数，/>分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)将图像重建的非线性问题转化为线性问题：边界测量电压与电导率分布的关系是非线性的，对于变化不大的电导率分布，可将边界测量电压的变化简化为线性形式式中，Δg是电导率分布的扰动，ΔU是电导率的变化引起的边界测量电压的变化，将以上线性化形式表示为b＝Ag；

(3)设置初始化参数：最小迭代阈值ε＝10^-6，初始迭代次数k＝0，最大迭代次数k_max＝200，惩罚参数δ₁＝200，δ₂＝δ₃＝30，增广拉格朗日乘子α⁰＝0，β⁰＝0，χ＝0初始灰度值g⁰＝A^Tb；

(4)更新权重因子ω：根据上一次优化得到的值自动选择合适的权重因子，

(5)更新非凸参数π₁，π₂：根据上次迭代结果值的一阶梯度和二阶梯度自动选择合适的非凸参数，/>

(6)更新正则化参数λ：根据上次迭代的值的最大值和最小值为先验经验自适应的选取合适的正则化参数，/>式中，r代表有限元剖分网格的平均长度，t为调节正则化参数的因子0＜t＜100；

(7)更新辅助变量v：

(8)更新辅助变量w：

(9)更新辅助变量z：

(10)更新电导率分布g：

为了提高计算效率和提高解的稳定性，采用快速傅里叶变换求解子问题：

式中，上标T是转置算子，和A^TA是具有循环块结构的块循环矩阵；

(11)更新增广拉格朗日乘子α,β,和χ：

(12)判断迭代是否符合迭代终止条件或者k＞k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1继续迭代求解；

(13)根据最终求解所得电导率值进行成像。