阅读说明：本技术 基于布尔运算的电阻抗成像方法 (Electrical impedance imaging method based on Boolean operation ) 是由 刘�东 顾丹萍 杜江峰 于 2020-11-11 设计创作，主要内容包括：一种基于布尔运算的电阻抗成像方法,包括以下步骤：建立完备电极模型,通过有限元方法求解得到观测模型；建立基于B样条曲线形式的基本单元形状的矩阵表示形式；对B样条曲线进行布尔运算；基于B样条曲线与布尔运算相结合的形状重构算法；执行重构算法,进行图像重构,以实现基于布尔运算的电阻抗成像。本发明所提出的方法可以重建得到多相、内嵌等复杂形状,得到的重建电阻抗图像成像精度、空间分辨率高,并且对电阻抗成像中普遍存在的模型误差、测量噪音等具有很好的鲁棒性。(An electrical impedance imaging method based on Boolean operation comprises the following steps: establishing a complete electrode model, and solving by a finite element method to obtain an observation model; establishing a matrix representation form of a basic unit shape based on a B spline curve form; performing Boolean operation on the B-spline curve; a shape reconstruction algorithm based on the combination of a B spline curve and Boolean operation; and executing a reconstruction algorithm to reconstruct the image so as to realize the electrical impedance imaging based on Boolean operation. The method provided by the invention can reconstruct and obtain complex shapes such as multiphase shapes, embedded shapes and the like, and the obtained reconstructed electrical impedance image has high imaging precision and spatial resolution and has good robustness to model errors, measurement noises and the like commonly existing in electrical impedance imaging.)

基于布尔运算的电阻抗成像方法

技术领域

本发明涉及电阻抗断层成像技术领域，是一种通过布尔运算与电阻抗成像技术相结合来重建对象体内阻抗分布图像的方法。

背景技术

电阻抗成像(Electrical impedance tomography，EIT)是一种无创、无辐射、低成本、能够连续监测的新型功能成像方法，此项医学影像技术是基于人体的各个组织具有不同电阻抗值这一特性，在人体外部施加一个比较小的安全激励电流(电压)，通过放置在体表的电极来测量电压(电流)，从而重建出人体内二维或三维的电阻抗分布图像。由于人体各个器官和组织的电导率不同，因此由EIT技术得到的电导率图像不仅包含了丰富的解剖学的信息，而且能反映出器官和组织电导率相应的生理、病理状态和功能信息。EIT技术在肺通气检测，肿瘤、癫痫等多种疾病的早期发现与诊断等方向有着诱人的应用前景。

自2019年12月起，世界人民面临着新型冠状肺炎的巨大考验，能够早期地发现疑似病人以及对确诊病人进行有效的治疗与监护对控制疫情的发展具有重要的意义。在《新型冠状病毒肺炎诊疗方案》(试行第七版)中提到严重感染者可快速进展为急性呼吸窘迫综合征，EIT可对感染者进行监护，并且在对病人采取肺复张治疗前可利用EIT技术进行可复张性评价。同时可通过EIT影像学监测肺炎感染者病情的发展，判断治疗的效果，最大限度地降低检查时往返于隔离室和放射科之间发生交叉感染的可能。EIT成像技术在新冠肺炎的监护、诊断等方面起着重要作用。

电阻抗技术具有无创、无辐射、功能信息丰富等优点，但其本质是一个病态的非线性、不适定的反问题，目前成像无法达到x射线计算机断层成像(CT)，核磁共振成像(MRI)等技术的成像精度，也由于EIT不适定的本质使得其对于多相、内嵌等复杂介质的重建能力有限，而在医学、工业应用中存在着大量复杂介质的重建情况。因此研究开发高精度、性能稳定的EIT系统及算法，提升其对于多相、内嵌等复杂介质的重建能力，探索其在临床医学以及工业领域中的应用，是当前的热点和难点问题。

在实现本发明的过程中，申请人发现了以下几个现有技术及其存在的缺点。

1、现有技术提出将B样条曲线应用于EIT。此算法将EIT图像重建问题转换为形状重建问题，利用B样条曲线通过调整控制点来逼近目标形状，提升了重建算法对于尖锐特征的保持能力，但是该算法并不具备拓扑演化能力，进行重建前需要获得区域内部物体的个数这一先验信息。对于无法知晓区域内部物体个数的情况，该算法则无法进行重建或无法重建得到可靠的重建结果，这在很大程度上限制这一算法的推广和应用。

2、现有技术中针对每一个目标物体给定多个MMC作为初始迭代值，通过优化组件的几何特征参数、改变组件之间的连通性或分布实现迭代过程形状的改变，并且利用高斯牛顿方法迭代改变电导率分布以及由正问题模型计算得到的电压值，从而通过逆问题最小化目标泛函，最终实现形状重建。但该发明仅针对两相的情况进行重建，对于多相以及内嵌等复杂形状无法进行重建。

3、现有技术中存在一种基于粒子群和正则化高斯-牛顿迭代算法的生物电阻抗成像方法。该方法在正问题计算中采用非均匀剖分，以提高成像精度。并且在生物电阻抗成像问题的求解中，采用粒子群算法产生一个接近真实值的初值作为正则化高斯-牛顿迭代算法的一个初始值，再利用正则化高斯-牛顿迭代算法求解逆问题。该方法采用非均匀剖分以及标准的粒子群法为正则化高斯-牛顿迭代算法产生一个初始值，可以提高成像的精确度并且克服牛顿类算法对初始值敏感的问题，虽然此方法充分利用了先验信息，但在一定程度上不仅增加了计算的复杂度，而且该方法重建的图像的空间分辨率，成像精度仍较低。

发明内容

有鉴于此，本发明的主要目的在于提供一种基于布尔运算的电阻抗成像方法，以期部分地解决上述技术问题中的至少之一。

为了实现上述目的，作为本发明的一方面，提供了一种基于布尔运算的电阻抗成像方法，包括以下步骤：

建立完备电极模型，通过有限元方法求解得到观测模型；

建立基于B样条曲线形式的基本单元形状的矩阵表示形式；

对B样条曲线进行布尔运算；

基于B样条曲线与布尔运算相结合的形状重构算法；

执行重构算法，进行图像重构，以实现基于布尔运算的电阻抗成像。

其中，所述方法适用于绝对成像、多相成像和差分成像。

其中，所述完备电极模型如下所示：

其中，σ(x)为电导率分布，x∈Ω为空间坐标，z_l为接触阻抗，U_l和I_l分别为电极e_l上的电压和电流，n表示外单位法向。

其中，所述观测模型的一般表达式为：

V＝U(σ)+e；

其中，V为测量电压，U(σ)为使用有限元法求解得到的正问题解，即计算电压，e为附加高斯噪音。

其中，所述建立基于B样条曲线形式的基本单元形状的矩阵表示形式包括：

将每一个基本单元边界均表示为B样条曲线形式；

将B样条曲线形式经过递归处理；

得到B样条曲线形式的矩阵表示形式。

其中，所述布尔运算对基本单元进行的操作包括联合、相交和相减。

其中，所述基于B样条曲线与布尔运算相结合的形状重构算法包括：

各基本单元经过布尔运算组成新的介质边界形状，而后对介质分别定义相应的电导率，电导率分布与B样条曲线的控制点相关，可表示为

其中，在所述电阻抗成像方法适用于绝对成像的情况下，所述观测模型具体表示为：

其中，V为测量电压，U为通过正问题求解所得的电压，e为附加噪声，δ为有限元网格离散程度参数，Ω为求解域。

其中，在所述电阻抗成像方法适用于绝对成像的情况下，基于最小二乘法及正则化技术，利用高斯牛顿方法求解下列表达式所示的最小化问题实现介质边界的形状重构：

式中，L_e为观测噪声协方差矩阵的Cholesky因子，满足Reg(·)为正则化项。

其中，所述高斯牛顿方法是通过迭代进行求解所述最小化问题的，迭代过程不断修正电导率值以及基本单元的控制点参数；和/或

在迭代过程中，需要求解雅克比矩阵

其中，能够由标准方法求解得到，能够通过扰动的方法进行计算；和/或

迭代过程采用线性搜索方法，直到迭代步长小于一个正值或者等于零时，迭代终止，得到最后的重建图像。

基于上述技术方案可知，本发明的基于布尔运算的电阻抗成像方法相对于现有技术至少具有如下有益效果之一或其中的一部分：

(1)利用B样条曲线表示基本单元，迭代过程中通过调整控制点的位置来对基本单元的形状进行局部的灵活调整，使算法具有重建物体细节特征的能力；

(2)将曲线的控制点作为设计变量，从而将图像重建问题转换为形状重建问题，有效降低未知数的数目以及计算时间与成本；

(3)结合布尔运算使算法具备拓扑演化能力，因此在进行重建之前无需知晓区域内所含的物体个数，从而使得算法能够适用于解决更多的实际问题；

(4)利用布尔运算对基本单元进行联合、相交、相减等操作可以构建得到多相、内嵌等复杂形状；

(5)本发明所提出的算法可以重建得到多相、内嵌等复杂形状，得到的重建电阻抗图像成像精度、空间分辨率高，并且对电阻抗成像中普遍存在的模型误差、测量噪音等具有很好的鲁棒性。

附图说明

图1是本发明实施例提供的联合、相减运算相结合的运算过程示意图；

图2是本发明实施例提供的方法流程图；

图3是本发明实施例提供的3次B样条闭合曲线示意图；

图4是本发明实施例提供的布尔运算后各区域电导率分布示意图；

图5是本发明实施例提供的算法流程图。

具体实施方式

随着科技的发展，人们生活水平不断地提高，人们开始更多地关注身体健康，越来越多的人希望得到更好的医疗服务。医学成像技术是各项医疗检查、医疗诊断中最重要的方式之一。目前常用的成像方式主要有MRI、CT、超声等，通过此类技术医生可以清楚的看到人体内部器官组织的成像，从而做出准确的诊断。但众所周知，CT、超声、MRI等也存在着不少的缺陷，不仅费用高昂，并且存在辐射，长期使用对人体有着非常大的伤害。

而电阻抗成像方法具有无创、无辐射、使用方便且设备价格相对低廉等优势，成为了近年来国内外研究的热点，在医学成像领域拥有着巨大的发展潜能，同时EIT成像技术在新冠肺炎的监护、诊断等方面起着重要作用。传统EIT算法，已经发展了30多年，但其成像精度、实时性、稳定性以及复杂形状重建等多方面仍存在不少局限性。而在医学和工业上大量存在重建区域内部复杂，对所含物体个数不明确的情况，因此提出一种具备拓扑演化能力，能够重建出多相、内嵌等复杂形状，且高分辨率和鲁棒性的电阻抗成像算法具有很大的实际需求。

本发明公开了一种基于布尔运算的电阻抗成像技术。该方法充分利用了B样条曲线的局部可修改性以及灵活控制的特性、布尔运算的拓扑演化以及复杂形状构造能力等优势，在原有的成像算法上进一步提升成像的质量，具备多相、内嵌等复杂形状的重建能力，并且对电阻抗成像中普遍存在的模型误差、测量噪音等具有很好的鲁棒性和稳定性。

具体地，本发明公开了一种基于布尔运算的电阻抗成像方法，包括以下步骤：

建立完备电极模型，通过有限元方法求解得到观测模型；

建立基于B样条曲线形式的基本单元形状的矩阵表示形式；

对B样条曲线进行布尔运算；

基于B样条曲线与布尔运算相结合的形状重构算法；

执行重构算法，进行图像重构，以实现基于布尔运算的电阻抗成像。

其中，所述方法适用于绝对成像、多相成像和差分成像。

其中，所述完备电极模型如下所示：

其中，σ(x)为电导率分布，x∈Ω为空间坐标，z_l为接触阻抗，U_l和I_l分别为电极e_l上的电压和电流，n表示外单位法向。

其中，所述观测模型的一般表达式为：

V＝U(σ)+e；

其中，V为测量电压，U(σ)为使用有限元法求解得到的正问题解，即计算电压，e为附加高斯噪音。

其中，所述建立基于B样条曲线形式的基本单元形状的矩阵表示形式包括：

将每一个基本单元边界均表示为B样条曲线形式；

将B样条曲线形式经过递归处理；

得到B样条曲线形式的矩阵表示形式。

其中，所述布尔运算对基本单元进行的操作包括联合、相交和相减。

其中，所述基于B样条曲线与布尔运算相结合的形状重构算法包括：

各基本单元经过布尔运算组成新的介质边界形状，而后对介质分别定义相应的电导率，电导率分布与B样条曲线的控制点相关，可表示为

其中，在所述电阻抗成像方法适用于绝对成像的情况下，所述观测模型具体表示为：

其中，V为测量电压，U为通过正问题求解所得的电压，e为附加噪声，δ为有限元网格离散程度参数，Ω为求解域。

其中，在所述电阻抗成像方法适用于绝对成像的情况下，基于最小二乘法及正则化技术，利用高斯牛顿方法求解下列表达式所示的最小化问题实现介质边界的形状重构：

式中，L_e为观测噪声协方差矩阵的Cholesky因子，满足Reg(·)为正则化项。

其中，所述高斯牛顿方法是通过迭代进行求解所述最小化问题的，迭代过程不断修正电导率值以及基本单元的控制点参数；和/或

在迭代过程中，需要求解雅克比矩阵

其中，能够由标准方法求解得到，能够通过扰动的方法进行计算；和/或

迭代过程采用线性搜索方法，直到迭代步长小于一个正值或者等于零时，迭代终止，得到最后的重建图像。

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明的目的是提出基于B样条曲线与布尔运算相结合的电阻抗成像方法，对医学、工业中存在的多相、内嵌等复杂介质进行重建，在一定程度上解决由于噪音和误差等所导致的电阻抗重建图像空间分辨率以及成像精度低的问题。

本发明给定多条B样条曲线作为初始迭代曲线，通过变动相应控制点的位置来调整B样条曲线所表示的基本单元的形状，并利用单个或多个组合形式的布尔运算来对基本单元进行操作，实现迭代过程图形形状的改变，利用高斯牛顿方法迭代改变电导率分布以及由正问题模型计算得到的电压值，从而通过逆问题求解最小化目标泛函，最终实现图像重建。

通过对基本单元进行联合、相交、相减等布尔运算操作构建得到复杂介质形状。以带有内嵌物质的肺部形状为例，可以通过多个形状联合运算组成肺部形状，再对其进行相减运算得到内嵌的物体，相关运算过程示意图如图1所示。

本发明提出基于布尔运算的电阻抗成像技术，不仅具有复杂形状的重建能力，而且对电阻抗成像中普遍存在的模型误差、测量噪音等均具有很好的鲁棒性。

本发明采用B样条曲线形式表示基本单元，并与布尔运算相结合来进行EIT形状重建，充分地利用了B样条曲线的局部可调整性以及布尔运算的拓扑演化等能力。该方法所具有的优势在于：①利用B样条曲线表示基本单元，通过调整控制点的位置来对基本单元进行局部的灵活控制，从而具有重建物体细节特征的能力；②将B样条曲线的控制点作为设计变量，很大程度上减少了未知数的维度以及计算成本和时间；③结合布尔运算使得算法具有了拓扑演化能力，重建时无需知晓区域内部物体的个数；④利用布尔运算对基本单元进行联合、相交、相减等操作可以构建得到多相、内嵌等复杂形状等。基于以上优势，产生了将其与电阻抗成像技术相结合的想法。

本发明主要针对EIT技术的逆问题进行了创新，初始给定多条B样条曲线(基本单元)，以此作为高斯牛顿迭代方法的初始值，通过求解目标泛函，对基本单元进行布尔运算操作(联合、相交、相减等)得到迭代形状，从而改变电导率分布以及计算电压值，最后实现形状重建。

本发明的流程示意图如图2所示。

本发明的具体技术方案如下：

1、正问题：对所需要求解区域进行有限元剖分，建立(FEM)有限元模型。

将16个电极均匀等距地放置在感兴趣区域边界上，轮流恒流激励其中的两个电极，引起区域内部的电压，通过建立完备电极模型(CEM)求解正问题，并在其余电极上测量相应的电压值。

CEM模型如下：

其中，σ(x)为电导率分布，x∈Ω为空间坐标，z_l为接触阻抗，U_l和I_l分别为电极e_l上的电压和电流，n表示外单位法向。

利用有限元法逼近可以得到正问题的观测模型：

V＝U(σ)+e

其中V为测量电压，U(σ)为使用有限元法求解得到的正问题解，即计算电压，e为加性高斯噪音，均值为e^*，协方差为Γ_e。

2、基于B样条曲线形式的基本单元形状隐式表达：

将每一个基本单元边界均表示为B样条曲线形式。以2维平面上的n+1个顶点作为B样条曲线的控制顶点，相应的均匀3次B样条曲线C可表示为B样条基函数的线性表示

其中，N_i，k(t)表示第i个k次(通常将k设为3)B样条基函数，通过递归的形式定义如下：

其中i＝0，1，…，n，并且表示均匀节点向量，定义为

B样条基函数与B样条曲线控制顶点是相互独立定义的，因此，B样条曲线可表示为矩阵相乘的形式

C＝NP

此处，N为B样条基函数所组成的矩阵，P为B样条曲线控制点横纵坐标所组成的矩阵。

如图3所示，为3次B样条闭合曲线示意图，图中，黑色的点为B样条曲线的控制点，黑色实线为B样条曲线，圆圈圈出的点表示重复控制点。

3、对B样条曲线进行布尔运算：

通过B样条曲线表达基本单元边界，定义水平集函数f(x)，在B样条曲线内f(x)＞0，在B样条曲线上f(x)＝0，在B样条曲线外f(x)＜0。由此，对于N_c条B样条曲线可得N_c个相对应的水平集函数{f_j，j＝1，2，…，N_c}。初始给定的B样条曲线可能是两两重叠，或是与其余B样条曲线均不重叠，而在优化过程当中，原先重叠的两条B样条曲线可能会不再重叠，原先不重叠的两条B样条曲线可能会变得相互重叠，这也是本发明所提出的算法的一大优势：拓扑演化能力。针对两两重叠的B样条曲线进行布尔运算相当于对水平集函数进行取max或取min计算。对多条B样条曲线可以采用单种布尔运算操作，同时也可以采用多种布尔运算组合形式对相应的水平集函数进行取max和取min的组合操作。

基于B样条曲线表示与布尔运算相结合的形状重构算法：

各基本单元经过布尔运算组成新的介质边界形状，而后对介质分别定义相应的电导率，下图给出了几种简要形状的电导率分布示意图。值得注意的是示意图中给出了单条B样条曲线情况、两条B样条曲线分别进行单种形式布尔运算以及多条B样条曲线进行多种组合形式布尔运算的情况。在实际算法设计当中，还存在着各式各样的形状以及布尔算法组合，这也体现了本发明所提出的电阻抗重建算法的灵活性。如图4所示，为布尔运算后各区域电导率分布示意图。

综上可知，电导率分布与B样条曲线的控制点相关，可表示为从数学理论上来看，电阻抗成像技术的求解属于求解椭圆型偏微分方程，在绝对成像框架下其观测模型可以表示为：

其中，V为测量电压，U为通过正问题求解所得的电压，e为附加噪声，δ为有限元网格离散程度参数，Ω为求解域。

在差分成像框架下其观测模型可以表示为：

ΔV≈JΔσ+Δe

其中，ΔV为前后两组测量电压的差值，为在初始背景电导率值假设下电压关于电导率的雅克比矩阵，为待求解的前后两组电导率差值，Δe为两组附加噪声之差。

为简便表示，此处以绝对成像为例，基于最小二乘法及正则化技术，利用高斯牛顿方法求解下列最小化问题实现介质边界的形状重构：

式中，L_e为观测噪声协方差矩阵的Cholesky因子，满足Reg(·)为正则化项。

通过高斯牛顿方法迭代求解上述最小化问题，迭代过程不断修正电导率值以及基本单元的控制点参数。

在迭代过程中，需要求解雅克比矩阵

其中，可由标准方法求解得到，可通过扰动的方法进行计算。

迭代过程采用线性搜索方法，直到迭代步长小于一个比较小的正值或者等于零时，迭代终止，得到最后的重建图像。如图5所示为算法的具体流程图。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。