CN111616708A - 一种精准识别脑卒中颅内病变区域的图像重建方法 - Google Patents

Abstract

本发明公开了一种精准识别脑卒中颅内病变区域的图像重建方法，用于电阻抗层析成像重建出血性脑卒中的电导率分布。该方法将电阻抗层析成像问题看作求解一个线性不适定问题，确定目标函数；重建方法的流程为：获取物场信息，求解相对边界测量值向量和灵敏度矩阵；设置初始化参数；使用加权矩阵转换数据保真项；利用迭代方式求解目标函数；判断迭代是否结束；根据求得的灰度值，进行成像。本发明有效的克服了边缘过光滑效应和抑制了阶梯伪影，提高了算法的稳定性、适用性和重建图像的空间分辨率，扩展了总变差正则化算法在医学检测领域的应用范围。

Description

一种精准识别脑卒中颅内病变区域的图像重建方法

技术领域

本发明属于电阻抗层析成像技术领域，具体涉及一种精准识别脑卒中颅内病变区域的图像重建方法。

背景技术

目前，出血性脑卒中被认为是一种严重的脑疾病，它会使脑组织暴露于血液的刺激作用下，导致患者的残疾或死亡。在出血性脑卒中的诊断中，计算机断层扫描(ComputedTomography,CT)和磁共振成像(Magnetic Resonance Imaging,MRI)是两种常用的影像学技术，具有较高的空间分辨率。然而，这两种技术都无法检测出血性中风的早期过程，阻碍了及时治疗。另外，CT和MRI的时间分辨率很低，成像时间较长，设备体积庞大，价格相对昂贵。

电阻抗层析成像(Electrical Impedance Tomography,EIT)技术是一种可以重建检测区域中的电导率分布技术，目前引起了人们的广泛关注。生物组织的生理或病理变化会引起电导率的变化，这可以通过用EIT获得的重建图像来反映。同时，与CT和MRI相比，EIT技术便携、安全、经济且无创。此外，还可以获得较高的时间分辨率。EIT技术由于其优越性，已被用于肺功能、乳腺癌、肝脂肪和脑功能的医学成像。特别是当出血性脑卒中发生时，由于血液的电导率高于正常脑组织的电导率，电导率的分布会发生变化，使得EIT甚至在出血性脑卒中的早期监测成为可能。然而，从EIT的测量数据恢复电导率分布通常是一个不适定反问题，不足以保证解的存在性、唯一性和稳定性，从而导致了重建图像的空间分辨率相对较低。

对于不适定问题，为了得到有效稳定的逆问题的解，通常选用正则化方法来逼近真实解。目前，典型的正则化方法有基于L2范数的Tikhonov正则化(TikhonovRegularization,TR)方法和基于L1范数的总变差(Total Variation，TV)正则化方法。Tikhonov正则化算法由于其简单性得到了广泛的应用，但图像重建质量依赖于正则化参数的选取，如M.Rezghi等人2009年发表于《计算和应用数学杂志》(Journal ofComputational and Applied Mathematics)第231卷，第914-924页，题为《基于Tikhonov正则化的L-曲线法的新的变形方法》(A new variant of L-curve for Tikhonovregularization)。Tikhonov正则化算法对于被测介质连续分布具有良好的性能。然而，当被测介质不连续分布时，在边界上施加了过度的光滑性，从而降低了重建图像的分辨率。为了保持锐利的边缘，总变差正则化算法提供了一种解决方案。该方法保留了边界的不连续性，并允许重建锐利的边缘，以产生更清晰的图像，在EIT领域受到了广泛的关注，如G González等人2017年发表于《计算机和数学及其应用》(Computers&Mathematics withApplications)第74卷，第564-576页，题为《电阻抗层析成像中的各向同性和各向异性总变差正则化》(Isotropic and anisotropic total variation regularization inelectrical impedance tomography)。然而，由于总变差方法往往会产生分段常数解，使得该方法在重建锐利边缘的同时也会在光滑区域产生“阶梯效应”，降低了重建图像的整体分辨率，严重限制了总变差方法在EIT中的应用。

为了减轻总变差正则化方法的阶梯效应，提高其重建图像的空间分辨率，本发明提出一种精准识别脑卒中颅内病变区域的图像重建方法来求解电阻抗层析成像的逆问题。通过比较该方法在被测区域内绝对重建和相对重建的图像效果，验证本发明所提算法的稳定性和适用性。

发明内容

本发明解决的技术问题是提供了一种具有较好稳定性和适用性的精准识别脑卒中颅内病变区域的图像重建方法，用于EIT重建出血性脑卒中的电导率分布。为了降低了求解目标函数的难度，提出一种加权矩阵对数据保真项进行转换；为了进一步去除伪影，提高鲁棒性，引入了软阈值算子；为了求解最小化的目标函数，采用了一种有效的迭代算法；为了验证所提方法重建图像的分辨率，使用数值模拟研究了绝对重建和相对重建的图像效果。与Tikhonov正则化方法和总变差正则化方法相比，本发明提出的精准识别脑卒中颅内病变区域的图像重建方法在电阻抗层析成像重建图像方面显示出了其有效抑制阶梯伪影和精准识别物场介质的优势。

本发明为解决上述技术问题采用如下技术方案，一种精准识别脑卒中颅内病变区域的图像重建方法，该方法将电阻抗层析成像作为一个线性不适定问题Ag＝b_meas，式中，A为灵敏度矩阵，b_meas为相对边界测量值向量，g为所求成像灰度值，所建立的最小化目标函数为：

其中

式中，S(g)是目标函数；λ是正则化参数，为一个正标量，用来保持保真度项和混合惩罚项之间的平衡；ω∈[0,1]是控制一阶TV和高阶TV惩罚项之间权衡的加权因子；||·||_p和||·||₁分别是向量Lp范数和L1范数；

表示梯度；l_x和l_y分别表示在x方向和y方向上的有限元剖分；

和

分别是一阶和二阶离散微分。

重建算法的流程如下：

本发明的有益效果是：本发明所述的精准识别脑卒中颅内病变区域的图像重建方法在计算目标函数时，为了保证算法的稀疏性和降低求解目标函数的难度，引入一种加权矩阵对目标函数进行转换；为了进一步去除伪影，提高鲁棒性，引入了软阈值算子；最后，为了求解最小化的目标函数，采用了一种有效地迭代算法。本发明提高了算法的稳定性和适用性，扩展了总变差正则化算法的应用范围，有效的克服了传统Tikhonov正则化算法的边缘过光滑效应和总变差正则化算法图像阶梯伪影严重的特点，在提高电阻抗层析成像技术重建图像的空间分辨率方面有很好的效果。

附图说明

图1为本发明的一种精准识别脑卒中颅内病变区域的图像重建方法的流程框图；

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

图3为本发明的有限元剖分网格(a)正问题有限元剖分(b)逆问题有限元剖分；

图4为本发明在选取六种模型(模型1-6)的真实分布时，典型的TR算法，TV算法和一种精准识别脑卒中颅内病变区域的图像重建方法，绝对重建的图像结果示意图；

图5为本发明在选取两组模型(模型7-12)的真实分布时，典型的TR算法，TV算法和一种精准识别脑卒中颅内病变区域的图像重建方法，相对重建的图像结果示意图；

图6为相对重建模型空间分辨率的定量评价。

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

具体实施方式

结合附图和实施例对本发明的一种精准识别脑卒中颅内病变区域的图像重建方法加以说明。

本发明的一种精准识别脑卒中颅内病变区域的图像重建方法，用于EIT重建出血性脑卒中的电导率分布。为了降低了求解目标函数的难度，提出一种加权矩阵对数据保真项进行转换；同时，为了进一步去除伪影，提高鲁棒性，引入了软阈值算子；最后，利用一种有效的迭代算法完成最终的逆问题求解。

如图1所示，为本发明的一种精准识别脑卒中颅内病变区域的图像重建方法流程图。

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

采取有限元剖分的方式来对EIT进行求解，如图3所示。其中，对于EIT正问题采用三角形网格剖分，如图3(a)所示；EIT逆问题采用正方形网格剖分，如图3(b)所示；为了避免图像重建过程中逆向违规对成像质量的影响，正问题的剖分网格数要大于逆问题的剖分网格数。

为了较好地体现本发明一种精准识别脑卒中颅内病变区域的图像重建方法与典型TR算法和TV正则化算法的不同，分别给出了六种绝对重建模型和两组相对重建模型在这三种正则化算法下的求解结果。

实施例包括如下具体步骤：一种精准识别脑卒中颅内病变区域的图像重建方法，该方法将电学层析成像作为一个线性不适定问题Ag＝b_meas。式中，A为灵敏度矩阵，b_meas为相对边界测量值向量，g为成像灰度值。

逆问题可以用最小二乘优化形式的目标函数表示：

式中，f(g)为目标函数。为了解决逆问题的不定性问题，正则化方法为提高解的稳定性提供了一种新的方法。在最小二乘的目标函数中增加正则化项，也称为惩罚项。正则化方法的一般形式可以描述为：

其中，λ是一个正则化因子，它控制最小二乘项(通常称为数据保真度项)

和正则化项R(g)之间的权衡。

TR算法可将正则化项R(g)替换为正则化项

可以描述为：

然而，TR算法在被测介质不连续分布时，在边界上施加了过度的光滑性，从而降低了重建图像的分辨率。为了保持锐利的边缘，TV正则化算法提供了一种解决方案。TV正则化算法的形式可以描述为：

虽然重建图像保持了边缘上不连续性，但会产生分段常数解，使得该方法在重建锐利边缘的同时也会在光滑区域产生“阶梯效应”，从而降低了重建图像的整体分辨率。

为了抑制TV正则化方法的阶梯效应和提高重建图像的分辨率，本发明针对数据保真项和正则化项，提出了一种精准识别脑卒中颅内病变区域的图像重建方法来求解电阻抗层析成像的逆问题。其最小化的目标函数为：

其中

式中，S(g)是目标函数；λ是正则化参数，为一个正标量，用来保持保真度项和混合惩罚项之间的平衡；ω∈[0,1]是控制一阶TV和高阶TV惩罚项之间权衡的加权因子；||·||_p和||·||₁分别是向量Lp范数和L1范数；

表示梯度；l_x和l_y分别表示在x方向和y方向上的有限元剖分；

和

分别是一阶和二阶离散微分。重建步骤如下：

(1)获取物场信息：

针对重建模型，分别获取各自重建所需的边界测量值和灵敏度矩阵。边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极(如图2所示)，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值；灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：

式中A_ij是第j个电极对对第i个电极对的灵敏度系数，φ_i，φ_j分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布。

(2)设置初始化参数：

正则化参数λ＝75，加权因子ω＝0.3，范数p＝0.9，正常数ε＝10^-6，迭代终止θ_min＝10^-4，最大迭代次数k_max＝50，惩罚参数β₁＝100，β₂＝6，β₃＝6，软阈值系数κ＝0.08，拉格朗日乘子

和初始解g⁰＝0^n×1。

(3)使用加权矩阵转换数据保真项：

为了降低了求解目标函数的难度，使用一种加权矩阵对数据保真项进行转换，则最小化目标函数可以转变为：

其中

式中，W^k是在上一次迭代中计算并在每次迭代中更新的加权矩阵，g^k+1和g^k是第(k+1)次解的列向量和第k次解的列向量，ε是用于防止分母为零的正小常数，选择为10^-6。

(4)对目标函数进行迭代求解：

通过引入多个辅助变量，将复杂约束问题分解为一系列简单约束子问题。同时，每个子问题都有闭式解，以便用现有的数值方法求解。为此，引入了三个辅助变量f、q和r，可表示：

其增广拉格朗日函数可以表示为：

式中，γ₁,γ₂和γ₃表示增广拉格朗日乘子，β₁,β₂,β₃>0是惩罚参数。因此，辅助变量和灰度值的解可以描述为：

对于f子问题可表示为：

式中，

是逐点积，“sgn”表示signum的函数；

对于q子问题可表示为：

对于r子问题可表示为：

对于灰度值g的子问题可以表示为：

转换可得：

式中，上标T表示转置运算符，

和A^TA为循环矩阵。

为了进一步减少阶梯效应，引入了一个感兴趣的区域，定义：

式中，R指感兴趣的区域，κ是区分背景和介质的阈值系数。因此，灰度值的解g可转化为：

g^k+1＝R.*g^k+1

式中，.表示的是元素与元素之间的对应相乘；

最后，拉格朗日乘子的计算和更新方式可以表示为：

(5)迭代终止准则：

判断迭代是否符合迭代终止条件

或者是否达到最大迭代次数k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步继续迭代求解。

(6)输出结果：

根据最终求解所得灰度值g，进行成像。

图4为本发明选取六种典型的介质模型，传统的TR算法，TV正则化算法和本方法，绝对重建图像的结果示意图。从重建图像中可以看出，用TR算法重建图像时，在重建图像的边界上观察到过大的平滑度，难以准确识别图像的大小和形状。与TR算法相比，TV方法在重建过程中避免了TR算法出现过多的平滑度，具有更好的重建效果。然而，在重建图像的背景中却发现了严重的阶梯效应。与TR和TV方法相比，本方法在重建过程中很好地消除了过度平滑和阶梯效应。重建图像的质量最高，背景清晰，边界清晰，能准确地检测出颅脑中出血组织的大小和位置。

在临床出血性脑卒中研究中，随着出血过程的发展，重建电导率分布的相对变化具有重要意义。相对重建的目的是获得两个时间间隔之间电导率分布情况的变化(模型9和模型12)，还描绘了模型7、模型8、模型10和模型11的绝对重建，以便清楚地了解算法的性能，如图5所示。很明显，TR方法在绝对重建中很难区分介质的边界，而在相对重建中则在重建图像的边界上观察到模糊现象。对于TV方法，在绝对重建和相对重建中，可以很好地识别边界，同时观察到严重的阶梯效应。另外，对于TR和TV两种方法，重建的介质都大于实际的变化。在绝对和相对重建图像中，采用本发明所提方法可以有效地抑制模糊，消除阶梯效应。此外，重构后的介质尺寸与实际尺寸几乎相同。同时，采用空间分辨率(SpatialResolution,SRE)来定量评价图像重建质量，其表达式可以表示为：

式中，

其中，其中P_d和P₀分别表示介质和整个检测区域的面积，P_i是第i个像素的面积，[g_d]_i＝1(如果g_i大于g的最大振幅的四分之一)，否则选择[g_d]_i＝0。

从图6可以看出，一种精准识别脑卒中颅内病变区域的图像重建方法与典型TR算法和TV正则化算法相比，具有最小的空间分辨率，能够有效的抑制边缘模糊，消除阶梯效应，并准确的描绘出被测场域内部介质分布，明显地提高了电阻抗层析成像逆问题求解精度和图像重建质量。

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

Claims (1)

1.一种精准识别脑卒中颅内病变区域的图像重建方法，其特征在于：将电阻抗层析成像作为一个线性不适定逆问题Ag＝b_meas，式中，A为灵敏度矩阵，b_meas为相对边界测量值向量，g为成像灰度值，建立最小化目标函数为：

其中

式中，S(g)是目标函数；λ是正则化参数，为一个正标量，用来保持保真度项和混合惩罚项之间的平衡；ω∈[0,1]是控制一阶TV和高阶TV惩罚项之间权衡的加权因子；||·||_p和||·||₁分别是向量Lp范数和L1范数；

表示梯度；l_x和l_y分别表示在x方向和y方向上的有限元剖分；

和

分别是一阶和二阶离散微分，具体重建步骤如下：

(1)获取物场信息：

针对重建模型分别获取重建所需的边界测量值和灵敏度矩阵，边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值；灵敏度矩阵是根据不含介质的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：

式中，A_ij是第j个电极对对第i个电极对的灵敏度系数，φ_i，φ_j分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：

正则化参数λ＝75，加权因子ω＝0.3，范数p＝0.9，正常数ε＝10^-6，迭代终止θ_min＝10^-4，最大迭代次数k_max＝50，惩罚参数β₁＝100，β₂＝6,β₃＝6，软阈值系数κ＝0.08，拉格朗日乘子

和初始解g⁰＝0^n×1；

(3)使用加权矩阵转换数据保真项：

为了降低求解目标函数的难度，使用一种加权矩阵对数据保真项进行转换，则最小化目标函数转变为：

其中，

式中，W^k是在上一次迭代中计算并在每次迭代中更新的加权矩阵，g^k+1和g^k是第(k+1)次解的列向量和第k次解的列向量，ε是用于防止分母为零的正小常数，选择为10^-6；

(4)对目标函数进行迭代求解：

通过引入多个辅助变量，将复杂约束问题分解为一系列简单约束子问题，同时每个子问题都有闭式解，以便用现有的数值方法求解，为此，引入了三个辅助变量f、q和r，具体表示为：

其增广拉格朗日函数表示为：

式中，γ₁,γ₂和γ₃是增广拉格朗日乘子，β₁,β₂,β₃>0是惩罚参数；

辅助变量和灰度值的解描述为：

对于f子问题表示为：

式中“°”是逐点积，“sgn”表示signum的函数；对于q子问题表示为：

对于r子问题表示为：

对于灰度值g的子问题表示为：

转换得到：

其中，上标T表示转置运算符，

和A^TA为循环矩阵；

为了进一步减少阶梯效应，引入一个感兴趣的区域，定义：

式中，R指感兴趣的区域，κ是区分背景和介质的阈值系数，灰度值的解g转化为：g^k+1＝R.*g^k+1，其中，.表示的是元素与元素之间的对应相乘；拉格朗日乘子的计算和更新方式表示为：

(5)迭代终止准则：

判断迭代是否符合迭代终止条件θ＝||g^k+1-g^k||₂/||g^k||₂≤θ_min或者是否达到最大迭代次数k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步继续迭代求解；

(6)输出结果：

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

Images