超声反射信息约束的腹部病变电阻抗图像重建方法与流程

本发明属于电阻抗层析成像技术领域，涉及实现一种采用超声反射提供的内含物边界信息作为约束的总变差正则化电阻抗层析成像方法并用于实现人体腹部器官内病灶的精确重建。

背景技术：

电阻抗层析成像技术(eit)是一种功能性成像技术，通过在被测场域表面布置电极阵列并施加一定的电流激励并得到边界电压数据，以此来重建被测场域内部的电导率分布情况。相比结构性成像技术例如计算机断层扫描成像(ct)和核磁共振成像(mri)，eit技术重建图像分辨率较低但成像速度大大加快，可以达到实时成像的要求。作为一种新型的医学影像技术，电阻抗层析成像具有无侵入、无辐射、体积小和成本低等优势，是一种理想的实时病情监测手段。与此同时，eit作为一种实时监测工具，在医学影像、流动监测、地质勘探、以及建筑结构检测等领域也有着广阔的发展前景。

eit成像分辨率较低的主要原因是其逆问题(即由边界测量值重建场内电导率分布)的过程具有严重的病态性，这意味着边界测量值的微小扰动会导致解的大幅度变化。同时，电学成像本身的非线性问题导致重建图像具有较多的伪影和噪声。为了克服eit逆问题的病态性和非线性，专家学者们提出了许多图像重建算法，其中的正则化方法是一种克服病态性的有效手段。该方法通过构造正则化惩罚项的方式将一定的先验信息融入到逆问题的求解过程中，约束解的搜索空间，指引解的优化方向，以达到改善病态性的目的。典型的正则化方法有m.vauhkonen等人1998年在《ieee医学成像期刊》(ieeetransactionsonmedicalimaging)第17卷，285-93页，发表的题为《在电阻抗成像中的tikhonov正则化和先验信息》(tikhonovregularizationandpriorinformationinelectricalimpedancetomography)的文章中提到的l2正则化方法、j.zhao等人在《ieee国际成像系统与技术会议》(ieeeinternationalconferenceonimagingsystemsandtechniques)第25-30页发表的题为《电阻层析成像中小物体成像的稀疏正则化方法》(sparseregularizationforsmallobjectsimagingwithelectricalresistancetomography)的文章中提到的l1正则化方法、a.borsic等人2007年在《逆问题》(inverseproblems)第99卷，a12-a12页发表的题为《电阻抗成像中的总变差正则化方法》(totalvariationregularizationinelectricalimpedancetomography)的文章中提出的总变差正则化(tv)方法等。不同的正则化项能够引入不同类型的先验信息，比如tikhonov先验的均匀分布信息、拉普拉斯先验的光滑性信息以及m.cheney等人在1990年在《国际成像系统与技术杂志》(internationaljournalofimagingsystems&technology)第2卷，第66-75页发表的《noser：一种求解逆电导率问题的算法》(noser:analgorithmforsolvingtheinverseconductivityproblem)中提出的noser先验对应的非均匀分布信息。

除病态性和非线性外，eit还存在场域中心敏感度低和重建内含物边界模糊两方面的缺陷。由于电场的分布特性和电学层析成像中激励电流的限制，eit技术的灵敏度分布在场域中心相对较低，不能有效对远离边界电极的内含物进行有效的重建。物体的形状往往是图像中的重要信息，其边界是一幅图像的重要特征，由于电学层析成像的软场特性，eit重建结果中内含物边界梯度低，电导率值平滑过渡，无法有效分辨出物体的轮廓。如d.liu等人2017年在《ieee医学成像汇刊》(ieeetransactionsonmedicalimaging)第9卷，第1页发表的文章《用于电阻抗成像的参数化水平集方法》(aparametriclevelsetmethodforelectricalimpedancetomography)，可以建立基于边界的电阻抗成像算法以重建内含物边界。然而在重建像素值的eit图像中，由于eit的低分辨率物体的边界往往相当模糊。利用内含物的形状先验信息，构造形状约束正则化项，有望能在eit重建图像中保留较明显的内含物边界。

人体组织结构复杂，不同组织和器官之间的电导率分布差别迥异，在对病灶进行术后监护或长期检测的过程中，需要对病灶尺寸和形状的改变有比较清晰的判断和了解。在医学成像领域中，形状先验信息的提供大多通过多模态成像方法融合来实现，如ali.f.m等人2010年在《现代光学杂志》(journalofmodernoptics)第57卷，第273-286页发表的《一种用于mr和ct图像融合的曲线转换方法》(acurvelettransformapproachforthefusionofmrandctimage)中提出的基于曲线波变换的ct-mri融合算法；guo等人2008年在《神经计算》(neurocomputing)第72卷，第203-211页发表的《基于多尺度几何分析和轮廓波变换的多模态医学图像融合》(multimodalitymedicalimagefusionbasedonmulti-scalegeometricanalysisofcontourlettransform)中提出了基于非采样轮廓波变换的mri-spect图像融合方法。同样作为一种无侵入、无辐射、体积小和成本低的成像方法，超声成像技术广泛应用于医学影像诊断和术后监护中，其使用环境同电学层析成像相似，特征相近；同时，超声层析成像技术具有硬场特性，对组织结构的变化相对比较敏感，其反射模态可以对病灶的形状和几何尺寸有较为清晰的判断，通过超声反射模态获取病灶边界信息并约束电阻抗成像的重建过程，可以将超声成像的结构成像优势和电学成像的功能成像优势结合起来，对被测病灶的尺寸及电导率变化均能够进行完整和清晰的表征。目前，超声/电学双模态方法大多采用各自模态分别成像并进行图像融合的策略，融合结果相对较差，不能充分利用超声的结构性先验信息。因此，需要一种基于超声反射的病灶边界约束信息来更准确的指导eit的重建过程。

技术实现要素：

本发明针对人体腹部病灶电阻抗层析成像图像重建中，传统电阻抗成像方法不能有效重建病灶尺寸及边界轮廓的问题，提出一种基于超声反射边界约束信息的电阻抗层析成像图像重建算法。

一种基于超声反射边界约束信息的电阻抗层析成像图像重建方法，包括如下步骤：

步骤一：根据被测场域，获取图像重建所需的边界电压测量差值δv，具体计算方式为

δv＝vmea-vref

式中vref表示通过仿真计算得到的参考场边界电压测量值，vmea为测得的存在内含物下的实际场边界电压测量值。

步骤二：根据电磁场的互易性质构建jacobian矩阵，基于邻域总变差正则化方法给出逆问题求解目标函数

[1]jacobian矩阵的获取，是指根据仿真计算得到的参考场边界电压测量值，结合互易定理理论，计算灵敏度矩阵，其计算公式为：

式中，sij表示第i个电极对相对第j个电极对的灵敏度系数，是jacobian矩阵s的第i行、第j列元素，φi，φj分别表示第i个电极对和第j个电极对在激励电流分别为ii和ij时的场域电势分布，表示梯度算子，∫x∫ydxdy表示对场域内各像素单元的长和宽进行积分。

[2]基于邻域总变差正则化方法给出电阻抗层析成像逆问题求解目标函数，其计算公式为：

其中，g表示重建图像结果中每个像素单元的电导率值，表示满足表达式取得最小值时g的取值，s表示jacobian矩阵，表示二范数的平方，λ表示总变差正则化的正则化参数，lp表示总变差正则化矩阵，通过对不同像素间位置关系计算得到，β表示一个事先选定的正常数，一般选取为0.01，其主要作用是防止当像素值梯度等于0时正则化项不可微的情况出现，p表示场域内第p个像素，重建图像像素单元总个数为n。

[3]将电阻抗层析成像逆问题求解目标函数使用最小二乘法进行展开，得到第k次迭代重建的目标函数，其计算公式为：

其中，g^k+1表示第k次迭代的目标函数，g^k表示第k次迭代采用的像素电导率值，δg^k表示第k次迭代所需优化的像素电导率值变化量。

步骤三：使用超声反射模态得到病灶边界轮廓点位置并构建约束方程：

[1]基于超声反射模态，超声换能器发射脉冲超声波，记录发射换能器及临近超声换能器上接收到的声压信号随时间的变化，记录发射声波和接收声波之间的渡越时间并计算边界轮廓点距探头或探头中心连线的垂直距离，计算方式可以写为：

d＝c·tf/2

其中d表示病灶边界轮廓点距探头或探头中心连线的垂直距离，c表示人体腹部组织的平均声速，tf表示发射声波和接收声波之间的渡越时间。根据边界轮廓点距探头或探头中心连线的垂直距离得到边界轮廓点的具体坐标。

[2]一般认为，腹部均质器官内病灶在二维界面下均为凸封闭曲线，为了定量化其形状指标，采用椭圆形状作为边界轮廓点的拟合目标，得到封闭的边界轮廓，表示为：

其中，x表示轮廓上点的横坐标，xc为拟合椭圆轮廓中心点的横坐标，y表示轮廓上点的纵坐标，yc为拟合椭圆轮廓中心点的纵坐标，a为拟合椭圆的长轴，b为拟合椭圆的短轴。

[3]构建边界点的约束方程，要求椭圆边界的梯度下降值取为较大的数值，其他区域的梯度下降值取为较小的数值。以此构建约束方程表示为：

h(g)＝[g(pi)-α·gb,...,g(pb)-gb,…,g(po)-α^-1·gb,…,g(p)-g]＝0

其中，pi表示拟合轮廓上像素内法向向量指向的像素单元，pb表示在拟合轮廓上的像素单元，po表示拟合轮廓上像素外法向向量指向的像素单元，gb表示拟合轮廓上像素单元的电导率值，其他的像素单元统一使用p表示，其电导率值使用g表示，α为梯度下降取值，选取为10⁵。

步骤四：基于拉格朗日乘子法对超声约束方程指导下的电阻抗层析成像目标函数进行优化求解。

[1]构建拉格朗日函数：在逆问题重建算法每次迭代过程中，电阻抗层析成像目标函数f(δg^k)可以表示为：

超声轮廓梯度约束方程g(δg^k)可以表示为：

g(δg^k)＝h(g^k)+jh(g^k)·δg^k＝0

其中jh(g^k)表示在第k次迭代过程中，方程h(g^k)的一阶偏微分矩阵。根据拉格朗日乘子法，构建新的拉格朗日函数l(δg^k,μ)，表示为：

l(δg^k,μ)＝f(δg^k)+μg(δg^k)

其中，μ为拉格朗日系数。

[2]拉格朗日函数求解：拉格朗日函数在取得极值(最小值)的条件必要条件为：

其中，表示偏微分求解。根据上述目标函数和约束方程，上式可以展开并表示为：

其中，

[3]根据拉格朗日乘子函数极值条件构建新的电学成像目标方程及优化求解目标函数，目标方程表示为：

将上式子使用snew·δx＝bnew表示，并使用高斯牛顿迭代对该方程进行求解，求解中每步迭代表示为：

δx^k+1＝δx^k-(snew^tsnew+ηi)^-1·snew^t·(snewδx^k-bnew)

其中，k表示迭代次数，i表示正则化矩阵(此处用单位阵替代)，η为高斯牛顿迭代中的正则化参数。

[4]通过上式高斯牛顿法的迭代，得到δg^k的权系数并用以更新各个像素单元的电导率分布，计算方式表示为：

g^k+1＝g^k+ξ·δg^k

其中，ξ为更新像素单元电导率值的步长。

步骤五：重复步骤四直至残差满足要求

其中，reg^k＝||s·g^k-δv||表示残差值,ε为人为设定的残差阈值。

本发明中所提出的方法，可以在电阻抗成像结果中保留较完整准确的内含物边界，在给出内含物位置准确位置、尺寸的同时有效减少了成像伪影，显著提高了腹部病灶eit图像逆问题求解的重建质量。

所提算法以邻域总变差正则化重建算法为基础，将超声反射得到的病灶轮廓边界点转变为像素梯度约束方程，指导图像重建算法目标函数的优化求解过程。所提方法解决了传统电学成像算法在进行人体腹部均质器官病灶可视化检测中边界模糊、伪影严重的问题，显著提高了电阻抗层析成像技术对不同尺寸、不同位置病灶的分辨能力。扩展了总变差正则化重建算法的应用，提高了电阻抗层析成像逆问题的求解精度和图像重建质量。所提出基于超声反射边界约束信息的电阻抗图像重建算法，其核心思想是“通过超声反射模态确定病灶边界轮廓点位置并转化为梯度约束方程，构造拉格朗日函数并求解重建电导率分布”其中：通过超声反射模态获取病灶边界轮廓点位置并拟合，实现了对病灶位置、尺寸等信息的有效获取；通过构建拉格朗日函数并使用高斯牛顿法求解，有效解决了传统电学成像图像重建算法对病灶重建结果边界模糊、尺寸分辨不清的问题。该算法可以在重建结果中保留较为清晰和准确的内含物形状结构，在保证成像速度的基础上显著提升电阻抗层析成像精度。

附图说明

图1为本发明的基于超声反射边界约束信息的电阻抗层析成像图像重建算法完整流程图，主要分为超声反射边界约束信息获取，与超声约束方程下电阻抗重建逆问题求解两部分。

图2为本发明中针对人体腹部病灶所使用的电阻抗层析成像系统结构图；

图3为本发明的五个典型仿真模型，并分别给出了相应的邻域总变差正则化(tv)成像结果及本发明算法的最终成像结果；

图4为针对本发明五组仿真模型不同成像结果的相对误差和图像相关系数对比。

具体实施方式

结合附图和实施例对本发明的基于超声反射边界约束信息的电阻抗层析成像图像重建算法加以说明。

本发明的基于超声反射边界约束信息的电阻抗层析成像图像重建算法，实施例中针对人体上腹部均质器官内病灶可视化监测这一可视化应用形式。使用人体上腹部二维截面结构先验构建表征肝脏、肾脏内中肿瘤电阻抗层析成像的正问题模型，通过超声反射模态确定病灶的边界轮廓位置信息并转化为约束方程。图像重建逆问题的迭代求解过程可以分解为基于邻域总变差正则化方法的拉格朗日函数构建和使用高斯牛顿法求解像素单元电导率迭代值两部分，拉格朗日函数的构建将超声约束方程和电学优化目标函数统一在统一的求解框架下，高斯牛顿法的迭代求解可以快速的在满足约束方程的条件下给出电导率迭代值的数值解，并用以更新像素单元的电导率分布和拉格朗日函数。

如图1所示，为本发明的基于超声反射边界约束信息的电阻抗层析成像图像重建算法。算法主要分为边界测量值获取、互易定理计算灵敏度、超声反射约束方程构建、拉格朗日函数构建，高斯牛顿法迭代计算电导率迭代值并更新拉格朗日函数五个部分。图2是其针对人体腹部内病灶的电阻抗层析成像测试方式的基本示意。

图3、图4中分别给出了传统电阻抗成像图像重建算法的成像结果与本算法的成像结果与重建指标对比，重建指标包括相对误差(re)和图像相关系数(cc)两种，其计算方法可以表示为：

其中，σ表示重建的像素单元电导率分布，σ^*表示真实情况下的电导率分布，σj和σj^*表示第j个像素单元重建的和真实的电导率分布，和表示重建的和真实的电导率分布的平均值。

本算法实施例包括如下具体步骤：

(1)针对图4中模型1-模型4的内含物分布情况，分别获取各自重建所需的边界电压测量数据，根据被测场域，获取图像重建所需的边界电压测量差值δv，具体计算方式为

δv＝vmea-vref

式中vref表示通过仿真计算得到的参考场边界电压测量值，vmea为测得的存在内含物下的实际场边界电压测量值，如图2所示，共32个电极参与测量与重建，根据“相邻电流激励，相邻电压测量”的激励采集策略，vmea、vref均包含928组电压测量数据。

(2)根据电磁场的互易性质构建jacobian矩阵，基于邻域总变差正则化方法给出逆问题求解目标函数

a.jacobian矩阵的获取，是指根据仿真计算得到的参考场边界电压测量值，结合互易定理理论，计算灵敏度矩阵，其计算公式为：

式中，sij表示第i个电极对相对第j个电极对的灵敏度系数，是jacobian矩阵s的第i行、第j列元素，φi，φj分别表示第i个电极对和第j个电极对在激励电流分别为ii和ij时的场域电势分布，表示梯度算子，∫x∫ydxdy表示对场域内各像素单元的长和宽进行积分。

b.基于邻域总变差正则化方法给出电阻抗层析成像逆问题求解目标函数，其计算公式为：

其中，g表示重建图像结果中每个像素单元的电导率值，表示满足表达式取得最小值时g的取值，s表示jacobian矩阵，表示二范数的平方，λ表示总变差正则化的正则化参数，选取经验值为3×10^-4，lp表示总变差正则化矩阵，通过对不同像素间位置关系计算得到，β表示一个事先选定的正常数,这里选取为0.01，其主要作用是防止当像素值梯度等于0时正则化项不可微的情况出现，p表示场域内第p个像素，重建图像像素单元总个数为n。

c.将电阻抗层析成像逆问题求解目标函数使用最小二乘法进行展开，得到第k次迭代重建的目标函数，其计算公式为：

其中，g^k+1表示第k次迭代的目标函数，g^k表示第k次迭代采用的像素电导率值，δg^k表示第k次迭代所需优化的像素电导率值变化量。

(3)使用超声反射模态得到病灶边界轮廓点位置并构建约束方程：

a.基于超声反射模态，超声换能器发射脉冲超声波，记录发射换能器及临近超声换能器上接收到的声压信号随时间的变化，记录发射声波和接收声波之间的渡越时间并计算边界轮廓点距探头或探头中心连线的垂直距离，计算方式表示为：

d＝c·tf/2

其中d表示病灶边界轮廓点距探头或探头中心连线的垂直距离，c表示人体腹部组织的平均声速，tf表示发射声波和接收声波之间的渡越时间。根据边界轮廓点距探头或探头中心连线的垂直距离得到边界轮廓点的具体坐标。

b.一般认为，腹部均质地器官内病灶在二维界面下均为凸封闭曲线，为了定量化其形状指标，采用椭圆形状作为边界轮廓点的拟合目标，得到封闭的边界轮廓，表示为：

其中，x表示轮廓上点的横坐标，xc为拟合椭圆轮廓中心点的横坐标，y表示轮廓上点的纵坐标，yc为拟合椭圆轮廓中心点的纵坐标，a为拟合椭圆的长轴，b为拟合椭圆的短轴。

c.构建边界点的约束方程，要求椭圆边界的梯度下降值取为较大的数值，其他区域的梯度下降值取为较小的数值。以此构建约束方程表示为：

h(g)＝[g(pi)-α·gb,...,g(pb)-gb,…,g(po)-α^-1·gb,…,g(p)-g]＝0

其中，pi表示拟合轮廓上像素内法向向量指向的像素单元，pb表示在拟合轮廓上的像素单元，po表示拟合轮廓上像素外法向向量指向的像素单元，gb表示拟合轮廓上像素单元的电导率值，其他的像素单元统一使用p表示，其电导率值使用g表示，α为梯度下降取值，选取为10⁵。

(4).基于拉格朗日乘子法对超声约束方程指导下的电阻抗层析成像目标函数进行优化求解。

a.构建拉格朗日函数：在逆问题重建算法每次迭代过程中，电阻抗层析成像目标函数f(δg^k)可以表示为：

超声轮廓梯度约束方程g(δg^k)可以表示为：

g(δg^k)＝h(g^k)+jh(g^k)·δg^k＝0

其中jh(g^k)表示在第k次迭代过程中，方程h(g^k)的一阶偏微分矩阵。根据拉格朗日乘子法，构建新的拉格朗日函数l(δg^k,μ)，表示为：

l(δg^k,μ)＝f(δg^k)+μg(δg^k)

其中，μ为拉格朗日系数,经验取值为5×10^-4。

b.拉格朗日函数求解：拉格朗日函数在取得极值(最小值)的条件必要条件为：

其中，表示偏微分求解。根据上述目标函数和约束方程，上式可以展开并表示为：

其中，

c.根据拉格朗日乘子函数极值条件构建新的电学成像目标方程及优化求解目标函数，目标方程表示为：

将上式子使用snew·δx＝bnew表示，并使用高斯牛顿迭代对该方程进行求解，求解中每步迭代表示为：

δx^k+1＝δx^k-(snew^tsnew+ηi)^-1·snew^t·(snewδx^k-bnew)

其中，k表示迭代次数，i表示正则化矩阵(此处用单位阵替代)，η为高斯牛顿迭代中的正则化参数，经验取值为3×10^-4。

d.通过上式高斯牛顿法的迭代，得到δg^k的权系数并用以更新各个像素单元的电导率分布，计算方式可以写为：

g^k+1＝g^k+ξ·δg^k

其中，ξ为更新像素单元电导率值的步长，经验取值为2.5×10^-2。

(5).更新拉格朗日函数并使用高斯牛顿迭代法求解，重复上述步骤直至残差满足要求：

其中，reg^k＝||s·g^k-δv||表示残差值,ε为人为设定的残差阈值，经验取值1×10^-4。

本发明中所提出的方法，可以在电阻抗成像结果中保留较完整准确的内含物边界，在给出内含物位置准确位置、病灶的同时有效减少了成像伪影，显著提高了腹部病灶eit图像逆问题求解的重建质量。通过超声反射模态获取病灶边界轮廓点位置并拟合，实现了对病灶位置、尺寸等信息的有效获取；通过构建拉格朗日函数并使用高斯牛顿法求解，有效解决了传统电学成像图像重建算法对病灶重建结果边界模糊、尺寸分辨不清的问题。该算法可以在重建结果中保留较为清晰和准确的内含物形状结构，在保证成像速度的基础上显著提升电阻抗层析成像精度。

以上所述实施例为本发明的几个示例模型，本发明不局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都在本发明保护的范围。