CN106037650A - 一种混合变差生物电阻抗成像方法 - Google Patents

Abstract

本发明一种混合变差生物电阻抗成像方法，涉及测量身体胸腔的电导率，该方法采用有限元方法来模拟电位空间分布情况，建立生物电阻抗成像技术的数值模型，将尺度较大型的数值模型离散为较小单元，施加一定的边界条件并计算边界电压，分别建立二维胸腔数值模型和三维胸腔数值模型，确定混合变差算法的目标函数后由最速下降法进行逆问题求解并进行胸腔电阻抗图像重建，以说明重建图像的分辨率的提高，克服了现有的生物电阻抗成像技术在重建分辨率更高的图像过程中受到噪声和误差的影响，出现过度光滑现象和阶梯效应，还不能获得准确的肺部电阻抗重建图像的缺陷。

Description

一种混合变差生物电阻抗成像方法

技术领域

本发明的技术方案涉及测量身体胸腔的电导率，具体地说是一种混合变差生物电阻抗成像方法。

背景技术

生物电阻抗成像技术(Electrical Impedance Tomography，EIT)是继形态、结构成像之后，于近三十年出现的新一代医学成像功能技术。根据人体内不同组织在不同生理、病理状态下具有不同的电阻抗的特点，通过体表电极给人体施加小的安全激励电流或电压，在人体表面测量响应的电压或电流信号，以重建人体内部的电阻抗分布图像。人体组织的生理功能变化或者病理改变会引起该组织电阻抗的变化，且这些变化能够在电阻抗成像的重构图像中显现出来，所以生物电阻抗成像是一种功能成像技术。由于生物电阻抗成像技术具有无创、无射线损伤、方便携带、系统结构简单、操作简便和适合于早期疾病(如肿瘤)筛查的特点，并以其时间分辨率高、成像速度快和功能成像的独特优势，已经初步应用于临床研究。生物电阻抗成像技术是当今生物电磁学与生物医学工程的重要研究课题之一。

然而，现有技术中应用生物电阻抗成像技术进行电特性成像时，所得图像的分辨率仍较低。这主要是因为成像过程中需要处理具有病态特征的逆问题。针对这一病态性，现有技术常采用吉洪诺夫算法进行求解。此类方法在对连续形式的电导率目标进行图像重建时，可以得到较好的重建效果；但是，在对具有跳跃式的电导率目标进行图像重建时，一般都会出现过光滑现象，会将跳跃部分的电导率光滑化，成像结果会存在伪影，丧失了生物电阻抗成像应有的高对比度和锐度，如器官与器官间边界上电特性的阶跃性变化，重建过程中往往会被吉洪诺夫算法平滑掉。总变差算法是一种可以不通过施加平滑进行正则化的方法，它能够消弱重构结果的连续性，可以保持重建图像中的非连续变化，增强间断性；但是，在处理逆问题过程中，如果目标函数具有连续性，采用总变差算法进行求解，则会出现阶梯效应，导致不能得到满意的重建效果。如人体肺部，既具有电导率变化明显的区域，也有电导率变化较缓和的区域，因此，如何获得准确的肺部重建图像以准确反映真实肺部电导率分布信息是目前亟需解决的问题。

通过对已有资料的检索，刘近贞的博士论文《基于扫描电极的开放式电阻抗成像的研究》中提到了混合变差算法，该论文中采用了牛顿-拉夫逊法对逆问题进行求解，而牛顿-拉夫逊方法的收敛性主要取决于一个良好的启动初值，如果初值选择不当，算法有可能不收敛，最终无法求出目标函数的最小值；该论文在选取正则化参数时，其选择依据为重建目标的大小和深度，由于目标的大小和深度是未知的，所以该方法有一定局限性。

另一方面，经过对现有生物电阻抗成像技术图像重建的检索，生物电阻抗成像技术成像实验都是在合理的生物电阻抗成像技术系统下进行的。由于生物电阻抗成像技术重建问题是严重病态的，且在图像重建过程中会受到仪器中误差和硬件中噪声源的影响，通过研究表明，仪器中的误差和硬件产生的噪声对重建的图像会产生较明显的影响，因此设计一套能减少噪声和误差影响的生物电阻抗成像技术系统对于重建出分辨率更高的图像是图像重建的关键步骤。

发明内容

本发明一种混合变差生物电阻抗成像方法，涉及测量身体胸腔的电导率，该方法采用有限元方法来模拟电位空间分布情况，建立生物电阻抗成像技术的数值模型，将尺度较大型的数值模型离散为较小单元，施加一定的边界条件并计算边界电压，分别建立二维胸腔数值模型和三维胸腔数值模型，确定混合变差算法的目标函数后由最速下降法进行逆问题求解并进行胸腔电阻抗图像重建，以说明重建图像的分辨率的提高，克服了现有的生物电阻抗成像技术在重建分辨率更高的图像过程中受到噪声和误差的影响，出现过度光滑现象和阶梯效应，还不能获得准确的肺部电阻抗重建图像的缺陷。

本发明解决该技术问题所采用的技术方案是：一种混合变差生物电阻抗成像方法，步骤如下：

A.设置进行混合变差生物电阻抗成像所用装置：

进行混合变差生物电阻抗成像所用装置的构成包括计算机与成像算法模块、通信模块、总控模块、信号发生模块、电荷泵电流源模块、激励通道选通模块、测量通道选通模块、电极阵列和信号调理模块；计算机与成像算法模块用来控制总控模块和程序运行，总控模块通过控制总线、数据总线和地址总线完成对生物电阻抗成像技术系统各部分精确实时控制和快速信息的交换，信号发生模块用于产生正弦电压信号，电荷泵电流源模块将电压转换成电流，激励通道选通模块用于选择和切换电极阵列中的不同电极注入模式以将电流信号注入被测目标，测量通道选通模块用于选择和切换电极阵列的测量模式，以提取被测目标表面感应的电压信号，送至信号调理模块，信号调理模块通过其内部的可编程增益放大电路、滤波电路、乘法解调电路和模数转换电路，对电压信号进行放大、滤波、解调和模数转换，转换为数字量的电压信号由通信模块送入总控模块，并由通信模块实时地送入计算机与成像算法模块，在计算机内调用成像算法完成图像的重建工作。

B.进行混合变差生物电阻抗成像方法以实现对胸腔的电阻抗图像重建，其技术方案如下：

第一步，正问题求解：

(1.1)建立基于有限元方法的生物电阻抗成像技术的数值模型(以下简称数值模型)；

(1.2)定义数值模型的物理特性；

(1.3)对数值模型进行离散；

(1.4)施加边界条件；

(1.5)计算边界电压；

第二步，确定算法的目标函数：

(2.1)确定吉洪诺夫算法的目标函数；

(2.2)确定总变差算法的目标函数；

(2.3)确定混合变差算法的目标函数；

第三步，进行逆问题求解和胸腔电阻抗图像重建：

首先进行逆问题求解，即采用最速下降法求解混合变差算法的目标函数；

然后，包括采用以下的胸腔电阻抗图像重建算法，

(3.1)采用吉洪诺夫算法进行胸腔电阻抗图像重建，

(3.2)采用总变差算法进行胸腔电阻抗图像重建，

(3.3)采用混合变差算法进行胸腔电阻抗图像重建，步骤如下：

(3.3.1)临床测得电压信号或实验仿体模型测得电压信号或由仿真计算得到电压数据；

(3.3.2)根据混合变差目标函数，应用优化的L曲线方法，自适应调节正则化参数λ²；

(3.3.3)重建胸腔电导率分布，获得重建的胸腔电阻抗图像；

第四步，输出图像：

由上述A所述的计算机输出重建胸腔电阻抗图像；

C.用上述A所述的进行混合变差生物电阻抗成像所用装置完成上述B所述的进行混合变差生物电阻抗成像方法以实现对胸腔的图像重建的技术方案的实施过程如下：

在计算机与成像算法模块输入确定用于激励人体胸腔的安全电流的幅值和频率，通过电极阵列将激励电流信号施加至人体胸腔表面对人体胸腔即被测目标进行激励，通过电极阵列测量胸腔表面感应的电压信号并送入计算机，在计算机内部，通过混合变差算法进行成像，总控模块为系统控制核心，接受来自计算机和成像算法模块的控制指令，实现对生物电阻抗成像技术硬件系统的整体控制和协调工作，信号发生模块产生1kHz-1MHz范围内可调的正弦波信号；总控模块控制整个系统的功能：在电荷泵电流源模块将电压信号转换成电流信号，电流信号的幅值位于0.1mA至1mA之间且可调节，通过激励通道选通模块，实现激励电极的选通与关闭，使电流信号按照设定的方式流入相应的电极，注入被测目标；通过控制测量通道选通模块选通电极阵列中的相应电极提取被测目标表面感应出的电压信号，以便将信号送入信号调理模块；控制信号调理模块内部的可编程增益放大电路的放大倍数，使得信号调理模块将感应电压信号进行有效放大后，再由信号调理模块内部的滤波电路和乘法解调电路分别对感应的电压信号进行滤波和解调，并由信号调理模块内部的模数转换电路将其转化为数字量信号；然后再将转化的数字量电压信号由通信模块送入计算机与成像算法模块，将数字量电压信号复原为模拟量电压信号，通过混合变差算法实现胸腔电阻抗图像重建，最后由计算机输出重建胸腔电阻抗图像。

上述一种混合变差生物电阻抗成像方法，所述根据混合变差目标函数，应用优化的L曲线方法，自适应调节正则化参数λ²，其具体的方法如下：

在优化的L曲线中，横坐标ξ是残差，纵坐标η是正则解，即

(ξ(λ²),η(λ²))＝(log||y-Ax||₂,log||x||₂)，

其中y为电压分布，A为雅克比矩阵，x为目标电导率分布；

对于混合变差正则化，优化的L曲线具有三个特征不同的部分，上部、中部和下部，其中上部和下部几乎垂直，中部几乎水平，优化的L曲线有两个拐角，其中，上部“垂直”部分对应着正则化参数过小时的解，下部“垂直”部分对应着正则化参数过大时的解，中部“水平”部分则对应着正则化参数取值适中时的解；优化的L曲线具有两个拐角，两个拐角将曲线的“水平”部分与上部“垂直”部分、下部“垂直”部分分开，通过定位两个拐角即可确定优化的L曲线的“水平”部分的范围，此范围对应着正则化误差和拟合误差之间达到平衡时的解；

上述定位拐角的方法是依据曲率最值准则，即选择优化的L曲线上曲率最大值点和最小值点，这两个拐角由下式求出：

$\mathrm{\κ}\left({\mathrm{\λ}}^2\right)=\arg \max \frac{{\mathrm{\ξ}}^{\′}{\mathrm{\η}}^{\′\′}-{\mathrm{\ξ}}^{\′\′}{\mathrm{\η}}^{\′}}{{({\left({\mathrm{\ξ}}^{\′}\right)}^2+{\left({\mathrm{\η}}^{\′}\right)}^2)}^{3/2}}$

其中ξ'、η'表示ξ、η的一阶导数，ξ″、η″表示ξ、η的二阶导数，

这样，由优化的L曲线“水平”部分来确定合适的正则化参数，

具体步骤如下：

(1)赋初值λ²＝10^-6，a＝0，计算混合变差算法目标函数的最小值，得到逆问题求解结果；

(2)将第一步的结果作为初值，进行下一步的迭代，根据第一步重构的结果，分别将a依次赋值为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1，针对每个给定的a值，依次寻找其对应的正则化参数λ²，并计算此时混合变差算法目标函数的最小值，重建电导率分布；当满足下式的停止准则时，迭代停止；

$\sqrt{\frac{\underset{i=1}{\overset{m}{\&Sigma;}}{\&lsqb;x\left(i\right)-x_0\left(i\right)\&rsqb;}^2}{\underset{i=1}{\overset{m}{\&Sigma;}}{\&lsqb;x\left(i\right)\&rsqb;}^2}}<{\mathrm{\&tau;\&delta;}}_n,$

其中，x(i)表示重建的电导率分布，x₀(i)表示设定的初始电导率分布，m表示电导率分布矢量的维数，τ为1.0至1.5之间的一个预设常数，δ_n为电压数据的噪声水平；

在利用最速下降法对混合变差算法的目标函数进行求解过程中，由于总变差算法的目标函数不是光滑函数，并非处处可导，本发明中为了解决▽x＝0时|▽x|不可微的情况，采用如下近似：

${||R\&CenterDot;\&dtri;x||}_1\&ap;{\&Integral;}_{\mathrm{\&Omega;}}\sqrt{{(R\&CenterDot;\&dtri;x)}^2+\mathrm{\&beta;}},$

其中R是正则化矩阵，Ω为求解场域，▽是梯度算子，β为一个选定的常数，则有：

$\begin{array}{c}{F}^{\&prime;}\left(x\right)=\{{(Ax-y)}^T(Ax-y)+{\mathrm{\&lambda;}}^2(1-a){\&lsqb;R(x-x_0)\&rsqb;}^T\&lsqb;R(x-x_0)\&rsqb;+{\mathrm{\&lambda;}}^2\&CenterDot;a{\&Integral;}_{\mathrm{\&Omega;}}\sqrt{{(R\&CenterDot;\&dtri;x)}^2+\mathrm{\&beta;}}d\mathrm{\&Omega;}\}\\ =x^T{A}^{T}A-y^{T}A+{\mathrm{\&lambda;}}^2\&lsqb;(1-a){(x-x_0)}^T\&CenterDot;R+a\&dtri;\&CenterDot;\frac{R\&CenterDot;\&dtri;x}{\sqrt{{\left|\&dtri;x\right|}^2}+\mathrm{\&beta;}}\&rsqb;\end{array},$

利用最速下降法求出混合变差算法的目标函数，求解中所依据的公式为：

x^k+1＝x^k-tF'(x)，

其中，t为最速下降法的步长，k是当前迭代的次数，每次迭代中采用线性搜索方法，使得不等式F^k+1(x)＜F^k(x)成立。

上述一种混合变差生物电阻抗成像方法，其中进行混合变差生物电阻抗成像所用装置是通过公知途径获得的，所涉及的操作方法是本领域技术人员所能掌握的。

本发明的有益效果是：与现有技术相比，本发明具有以下突出的实质性特点和显著进步：

(1)针对刘近贞的博士论文在选取正则化参数时，其选择依据为重建目标的大小和深度，由于目标的大小和深度是未知的，所以该方法有一定局限性，为了突破这一局限，本发明采用优化的L曲线方法自适应地调节正则化参数，同时，为了确保可以求出目标函数的最小值，本发明采用了初始值对收敛性影响较低的最速下降法来求解目标函数的最小值。

(2)本发明方法采用有限元方法来模拟电位空间分布情况，建立生物电阻抗成像技术的数值模型，将尺度较大型的数值模型离散为较小的对象，施加一定的边界条件并计算边界电压，分别建立二维胸腔数值模型和三维胸腔数值模型，通过确定混合变差算法的目标函数来进行逆问题求解并进行胸腔电阻抗图像重建，以说明重建图像的分辨率的提高，克服了现有的生物电阻抗成像技术在重建分辨率更高的图像过程中受到噪声和误差的影响，还不能获得准确的肺部电阻抗重建图像的缺陷。

(3)本发明在对比研究了现有的吉洪诺夫算法和总变差算法的基础上，提出同时具备这两种算法优势的混合变差算法，采用将吉洪诺夫算法和总变差算法融合在一起的混合变差算法对逆问题进行求解，该算法包括两部分，第一部分为拟合项，第二部分为正则化项，其中，正则化项将吉洪诺夫算法和总变差算法融合在一起，减少了计算量、提高了求解速度、增强了收敛稳定性，还保存了肺部电导率尖锐变化的部分，减少了由数据异常值带来的误差。这种算法不仅保留了吉洪诺夫算法在重建中的稳定性和良好的收敛性，同时也具备总变差算法识别目标边界的能力，提高了图像的分辨率，有效地改善了图像质量，可以得到较好的胸腔电阻抗重建图像。稳定且高分辨率的图像重建，也有助于一些与疾病(如肿瘤)相关的生物组织电特性的检测。可见，本发明的混合变差算法在胸腔电阻抗图像重建与临床研究中具有很重要的意义。

(4)本发明方法采用优化的L曲线方法确定合适的参数，其中，拐角将曲线的“水平”部分与两个“垂直”部分分开，“水平”部分对应着正则化误差和拟合误差之间达到平衡时的解，此范围内相应的参数是合适的正则化参数。在含有噪声、异常值等因素干扰的情况下，为逆问题求解提供了灵活性。

(5)本发明方法采用最速下降法来求解目标函数的最小值，该方法只需求一阶导数，所占的存储单元少，初始值对收敛性影响很低且收敛速度较快。

(6)本发明方法通过混合变差算法重建的图像的分辨率比采用吉洪诺夫算法重建的图像的分辨率提高了4.23％,比采用总变差算法重建的图像的分辨率提高了5.65％，对于本发明中的二维胸腔数值模型算例，重建速度比采用吉洪诺夫算法的重建速度提高了8.68％，比采用总变差算法的重建速度提高了10.91％。

(7)本发明方法的混合变差方法，不仅可以解决二维胸腔数值模型电导率重建问题，还可以解决具有严重病态特征的三维胸腔数值模型电导率重建问题；本发明既可用于临床实测数据的图像重建，也可用于仿真研究与实验中；本发明还可以对变形后的目标进行图像重建。

(8)本发明方法在建立的胸腔生物电阻抗成像模型的基础上进行图像重建，使得结果更加清晰，更有利于分析应用混合变差算法进行肺部电阻抗图像重建的优点。

(9)本发明方法也可应用于岩土工程、资源和环境保护领域中检测对象，该成像方法尤其适用于针对电特性变化进行的检测与成像的场合。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1是本发明方法实现对胸腔电阻抗图像重建的技术方案的总流程示意图。

图1.1是图1所示总流程中的正问题求解的流程示意图。

图1.2是图1所示总流程中的确定算法目标函数的流程示意图。

图1.3是图1所示总流程中的进行逆问题求解和胸腔电阻抗图像重建的流程示意图。

图2是本发明的进行逆问题求解和胸腔电阻抗图像重建所应用优化的L曲线示意图。

图3是现有技术中采用有限元方法建立的胸腔数值模型的图像。

图4是现有技术中吉洪诺夫算法重建的胸腔数值模型的图像。

图5是现有技术中总变差算法重建的胸腔数值模型的图像。

图6是本发明的混合变差算法重建的胸腔数值模型的图像。

图7是本发明的混合变差算法建立的三维胸腔数值模型的图像，其中：

图7(a)为混合变差算法建立的三维胸腔数值模型的图像的70度角正视图；

图7(b)为混合变差算法建立的三维胸腔数值模型的图像的90度角俯视图。

图8是本发明方法进行混合变差生物电阻抗成像所用装置的结构示意图。

图9是本发明方法完成的人体肺部呼吸电阻抗图像重建的序列图像。

图中，1.计算机与成像算法模块，2.通信模块，3.总控模块，4.信号发生模块，5.电荷泵电流源模块，6.激励通道选通模块，7.测量通道选通模块，8.电极阵列，9.信号调理模块,10.被测目标。

具体实施方式

图1所示实施例表明，本发明方法实现对胸腔电阻抗图像重建的技术方案的总流程是：正问题求解→确定算法目标函数→进行逆问题求解和胸腔电阻抗图像重建→输出图像→结束。

图1.1所示实施例表明，图1所示总流程中的正问题求解的流程是：建立基于有限元方法的生物电阻抗成像技术的数值模型→定义数值模型的物理特性→对数值模型进行离散→施加边界条件→计算边界电压。

图1.2所示实施例表明，图1所示总流程中的确定算法目标函数的流程是：确定吉洪诺夫算法的目标函数→确定总变差算法的目标函数→确定混合变差算法的目标函数。

图1.3所示实施例表明，图1所示总流程中的进行逆问题求解和胸腔电阻抗图像重建的流程是：首先进行逆问题求解，即采用最速下降法求解混合变差算法的目标函数→采用以下的胸腔电阻抗图像重建算法：采用吉洪诺夫算法进行胸腔电阻抗图像重建→采用总变差算法进行胸腔电阻抗图像重建→采用混合变差算法进行胸腔电阻抗图像重建。

图2显示本发明的进行逆问题求解和胸腔电阻抗图像重建所采用优化的L曲线，该曲线表明，本发明采用的优化的L曲线被两个拐角分成三部分，即上部“垂直”部分、“水平”部分和下部“垂直”部分，曲线特征明显，有利于本发明混合变差正则化参数合理范围的确定。

图3显示现有技术中采用有限元法对于任意形状物体建模具有很大优势，本发明中采用有限元方法建立的胸腔数值模型，该模型用于本发明正问题求解中计算边界电压数据和本发明逆问题求解与胸腔电阻抗图像重建。

图4显示现有技术中吉洪诺夫算法重建的胸腔数值模型的图像，由于吉洪诺夫算法对边界的过度光滑性，使得重建的胸腔电阻抗图像具有较大的图像伪影。

图5显示现有技术中总变差算法重建的胸腔数值模型的图像，该图像的边界具有明显的阶梯效应，不能得到准确的胸腔阻抗重建图像。

图6显示本发明的混合变差算法重建的胸腔数值模型的图像，该图像表明本发明采用的混合变差算法不仅克服了吉洪诺夫算法的过光滑性，也避免了总变差算法的阶梯效应，可以得到分辨率较高的胸腔电阻抗重建图像。

图7显示本发明的混合变差算法建立的三维胸腔数值模型的图像，该图像表明本发明建立的三维胸腔数值模型的图像比二维胸腔圆模型更精确，这证明本发明所采用的混合变差算法可以提高胸腔电阻抗重建图像的分辨率，更有利用分析胸腔呼吸时阻抗变化的胸腔电阻抗重建图像。

图8所示实施例表明，本发明方法进行混合变差生物电阻抗成像所用装置的结构是：计算机与成像算法模块1、通信模块2、总控模块3、信号发生模块4、电荷泵电流源模块5、激励通道选通模块6、测量通道选通模块7、电极阵列8和信号调理模块9；计算机与成像算法模块1用来控制总控模块3和程序运行，总控模块3通过控制总线、数据总线和地址总线完成对生物电阻抗成像技术系统各部分精确实时控制和快速信息的交换，信号发生模块4用于产生正弦曲线的电压信号，电荷泵电流源模块5将电压转换成电流，激励通道选通模块6用于选择和切换电极阵列8中的不同电极注入模式将电流信号注入被测目标10，测量通道选通模块7用于选择和切换电极阵列8的测量模式，电极阵列8位于被测目标的表面，以提取被测目标10表面感应的电压信号，送至信号调理模块9，信号调理模块9通过其内部的可编程增益放大电路、滤波电路、乘法解调电路和模数转换电路，对电压信号进行放大、滤波、解调和模数转换，转换为数字量的电压信号由通信模块2送入总控模块3，并由通信模块2实时地送入计算机与成像算法模块1，在计算机内调用成像算法完成图像的重建工作。

图8所示实施例还表明，本发明方法进行混合变差生物电阻抗成像所用装置的硬件部分和软件是相互兼容的，可以通过该装置测得人体胸腔的实际数据并进行图像重建，获得人体真实胸腔的图像重建结果。

图9显示本发明方法完成的人体肺部呼吸电阻抗图像重建的序列图像，该图像清晰显示呼吸过程中肺部阻抗分布的实时变化情况，表明本发明可以完成人体胸腔电阻抗图像的实时重构和显示。

实施例

本实施例在进行二维正问题仿真实验中，采用归一化的二维圆域模型，该模型模拟人的胸腔结构，在注入电流时求解胸腔电位的分布情况，其中该圆的半径为15cm。将归一化的圆域模型剖分成比较均匀的三角形网格，电极点位于圆边缘上，电极个数为16个，用16个圆点表示，注入50kHz，1mA的激励电流，设模型具有各向同性电导率分布，背景电导率为1Ω·m，目标电导率设为2Ω·m。本实施例采用有限元方法建立的胸腔二维正问题模型，其中单元数总数为6400个，节点总数为3281个。

本实施例在进行三维正问题仿真实验中，采用胸腔模型，模型采用64个圆形电极，电极置于胸腔模型的侧表面，该模型的高度为30cm，其中电极既可用于施加电流激励，建立三维敏感场，又能测量模型表面的电压信号，对其进行有限元剖分，剖分为一系列的四面体单元，设模型具有各向同性的电导率分布，背景电导率的值为1s/m，目标电导率为2s/m，注入50kHz,0.2mA的激励电流。本实施例采用有限元法建立的胸腔三维正问题模型，其中单元总数为126192个，节点总数为25731个,剖分层数为18层。

本实施例的一种混合变差生物电阻抗成像方法，步骤如下：

A.设置进行混合变差生物电阻抗成像所用装置：

本实施例进行混合变差生物电阻抗成像所用装置如图7所示实施例；

B.进行混合变差生物电阻抗成像方法以实现对胸腔的图像重建，其技术方案如下：

第一步，正问题求解：

(1.1)建立基于有限元方法的生物电阻抗成像技术的数值模型(以下简称数值模型)：

首先将电阻抗成像问题近似为稳态电流场来处理，等效为如下的Laplace方程的边值问题，

其中，σ为目标内部的电导率分布，为节点电位，

其次建立近似数值模型，其中包括具有规则形状的圆形数值模型，然后建立二维胸腔数值模型，最后建立三维胸腔数值模型；

(1.2)定义数值模型的物理特性：

定义的物理特性包括模型的尺寸、介质的电导率、电极的尺寸、电极个数、电极形状、电极的接触阻抗、注入电流的频率和幅值；

二维胸腔数值模型中的物理特性包括：1)采用归一化的二维圆模型，该模型模拟人的胸腔结构；2)该圆的半径为15cm；3)点状电极位于圆模型的边界节点上，电极个数为16个，电极的半径为0.15cm；4)注入50kHz，1mA的正弦波形激励电流；5)每个电极的接触阻抗为50Ω；6)设模型背景具有各向同性的电导率分布，电导率的值为1Ω·m,成像目标区域电导率为2Ω·m；

三维胸腔数值模型中的物理特性，1)利用胸腔CT图像勾勒出胸腔和肺部轮廓，构建真实尺寸的胸腔数值模型；2)模型的高度为30cm；3)模型采用64个圆形电极，电极分4层布置，每层16个电极，每层电极电极均匀布置于胸腔数值模型的表面，电极的半径为0.15cm；4)每个电极的接触阻抗为50Ω；5)每层电极按照相邻激励-相邻测量的模式，4层同时注入50kHz，0.2mA的正弦波形激励电流；6)设模型背景具有各向同性的电导率分布，电导率的值为1Ω·m,肺部所在区域电导率为2Ω·m；

(1.3)对数值模型进行离散：

利用有限元方法将尺度相对较大的数值模型离散为较小的单元，离散后二维胸腔数值模型单元总数为6400个，三维胸腔数值模型单元总数为126192个，并对各个小单元对象进行编号、记录并存储其节点所在位置；

(1.4)施加边界条件：

边界条件为

其中为边界上的电位，J_n为边界上注入电流的电流密度；

(1.5)计算边界电压：

由Laplace方程相应的泛函完成边界电压的计算，所用的Laplace方程泛函为

对泛函进行求解，求得电位的值，▽是梯度算子；

第二步，确定算法的目标函数：

(2.1)确定吉洪诺夫算法的目标函数：

首先，将电阻抗成像近似为如下线性系统

Ax＝y

其中，x为目标电导率分布，A为雅克比矩阵，y为电压分布；

常规的解决办法就是线性最小二乘法，即求解

$F\left(x\right)=\arg \underset{x}{min}{||Ax-y||}_2^2,$

其中||·||是欧几里得范数。该欠定问题的解不唯一，为了获得稳定的解，将最小二乘问题转化为

$F\left(x\right)=\arg \underset{x}{min}\{{||Ax-y||}_2^2+{\mathrm{\&lambda;}}^2{||R(x-x_0)||}_{L_P}^2\},$

其中，为拟合项，为正则化项，λ²是正则化参数，x₀是平均电导率的先验估计，当P的取值为2时，该函数是吉洪诺夫算法的目标函数；

(2.2)确定总变差算法的目标函数：

当P的值取1时，函数F(x)为总变差算法的目标函数；

(2.3)确定混合变差算法的目标函数：

为了避免总变差算法的阶梯现象与吉洪诺夫算法过度光滑效应，将求解最小化问题转化为

$F\left(x\right)=\arg \underset{x}{min}\{{||y-Ax||}_2^2+{\mathrm{\&lambda;}}^2(a{||R\&dtri;x||}_1+(1-a\left){||R\left(x-x_0\right)||}_2^2\right)\},$

其中，为拟合项，为混合变差正则化项，其中λ²(a||R▽x||₁)为总变差正则化分量，为吉洪诺夫正则化分量；a是权重函数，a的取值位于0到1之间；▽x＝x-x₀，吉洪诺夫算法和总变差算法之间的权衡关系是通过正则化参数λ²和a来确定；

第三步，进行逆问题求解和胸腔电阻抗图像重建：

首先进行逆问题求解，即采用最速下降法求解混合变差算法的目标函数；

然后，包括采用以下的胸腔电阻抗图像重建算法，

(3.1)采用吉洪诺夫算法进行胸腔电阻抗图像重建，

(3.2)采用总变差算法进行胸腔电阻抗图像重建，

(3.3)采用混合变差算法进行胸腔电阻抗图像重建，步骤如下：

(3.3.1)临床测得电压信号或实验仿体模型测得电压信号或由仿真计算得到电压数据；

(3.3.2)根据混合变差目标函数，应用优化的L曲线方法，自适应调节正则化参数λ²，

所述根据混合变差目标函数，应用优化的L曲线方法，自适应调节正则化参数λ²，具体的方法如下：

在优化的L曲线中，横坐标ξ是残差，纵坐标η是正则解，即

(ξ(λ²),η(λ²))＝(log||y-Ax||₂,log||x||₂)，

其中，y为电压分布，A为雅克比矩阵，x为目标电导率分布；

对于混合变差正则化，优化的L曲线具有三个特征不同的部分，上部、中部和下部，其中上部和下部几乎垂直，中部几乎水平，优化的L曲线有两个拐角，其中，上部“垂直”部分对应着正则化参数过小时的解，下部“垂直”部分对应着正则化参数过大时的解，中部“水平”部分则对应着正则化参数取值适中时的解；优化的L曲线具有两个拐角，两个拐角将曲线的“水平”部分与上部“垂直”部分、下部“垂直”部分分开，通过定位两个拐角即可确定优化的L曲线的“水平”部分的范围，此范围对应着正则化误差和拟合误差之间达到平衡时的解；

上述定位拐角的方法是依据曲率最值准则，即选择优化的L曲线上曲率最大值点和最小值点，这两个拐角由下式求出：

$\mathrm{\&kappa;}\left({\mathrm{\&lambda;}}^2\right)=\arg \max \frac{{\mathrm{\&xi;}}^{\&prime;}{\mathrm{\&eta;}}^{\&prime;\&prime;}-{\mathrm{\&xi;}}^{\&prime;\&prime;}{\mathrm{\&eta;}}^{\&prime;}}{{({\left({\mathrm{\&xi;}}^{\&prime;}\right)}^2+{\left({\mathrm{\&eta;}}^{\&prime;}\right)}^2)}^{3/2}}$

其中,ξ'、η'表示ξ、η的一阶导数，ξ″、η″表示ξ、η的二阶导数，

这样，由优化的L曲线“水平”部分来确定合适的正则化参数，

具体步骤如下：

(1)赋初值λ²＝10^-6，a＝0，计算混合变差算法目标函数的最小值，得到逆问题求解结果；

(2)将第一步的结果作为初值，进行下一步的迭代，根据第一步重构的结果，分别将a依次赋值为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1，针对每个给定的a值，依次寻找其对应的正则化参数λ²，并计算此时混合变差算法目标函数的最小值，重建电导率分布；当满足下式的停止准则时，迭代停止；

$\sqrt{\frac{\underset{i=1}{\overset{m}{\&Sigma;}}{\&lsqb;x\left(i\right)-x_0\left(i\right)\&rsqb;}^2}{\underset{i=1}{\overset{m}{\&Sigma;}}{\&lsqb;x\left(i\right)\&rsqb;}^2}}<{\mathrm{\&tau;\&delta;}}_n,$

其中，x(i)表示重建的电导率分布，x₀(i)表示设定的初始电导率分布，m表示电导率分布矢量的维数，τ为1.0至1.5之间的一个预设常数，δ_n为电压数据的噪声水平；

在利用最速下降法对混合变差算法的目标函数进行求解过程中，由于总变差算法的目标函数不是光滑函数，并非处处可导，本发明中为了解决▽x＝0时|▽x|不可微的情况，采用如下近似：

${||R\&CenterDot;\&dtri;x||}_1\&ap;{\&Integral;}_{\mathrm{\&Omega;}}\sqrt{{(R\&CenterDot;\&dtri;x)}^2+\mathrm{\&beta;}},$

其中R是正则化矩阵，Ω为求解场域，▽是梯度算子，β为一个选定的常数，则有：

$\begin{array}{c}{F}^{\&prime;}\left(x\right)=\{{(Ax-y)}^T(Ax-y)+{\mathrm{\&lambda;}}^2(1-a){\&lsqb;R(x-x_0)\&rsqb;}^T\&lsqb;R(x-x_0)\&rsqb;+{\mathrm{\&lambda;}}^2\&CenterDot;a{\&Integral;}_{\mathrm{\&Omega;}}\sqrt{{(R\&CenterDot;\&dtri;x)}^2+\mathrm{\&beta;}}d\mathrm{\&Omega;}\}\\ =x^T{A}^{T}A-y^{T}A+{\mathrm{\&lambda;}}^2\&lsqb;(1-a){(x-x_0)}^T\&CenterDot;R+a\&dtri;\&CenterDot;\frac{R\&CenterDot;\&dtri;x}{\sqrt{{\left|\&dtri;x\right|}^2}+\mathrm{\&beta;}}\&rsqb;\end{array},$

利用最速下降法求出混合变差算法的目标函数，求解中所依据的公式为：

x^k+1＝x^k-tF'(x)，

其中，t为最速下降法的步长，k是当前迭代的次数，每次迭代中采用线性搜索方法，使得不等式F^k+1(x)＜F^k(x)成立；

(3.3.3)获取胸腔电导率分布，输出重建的胸腔电阻抗图像；

第四步，输出图像：

由上述A所述的计算机输出重建胸腔电阻抗图像。

C.用上述A所述的进行混合变差生物电阻抗成像所用装置完成上述B所述的进行混合变差生物电阻抗成像方法以实现对胸腔的图像重建的技术方案的实施过程如下：

在人体体表剑突骨下沿以下2cm至剑突骨下沿以上13cm共计15cm范围内上下等间距均匀布置4层ECG电极，每层16个，共64个，采用本实施例进行混合变差生物电阻抗成像所用装置获取人体胸腔的实测数据，根据测得的实际数据进行图像重建；本实验分别获取了5组成年男性胸腔的真实数据，并利用混合变差算法进行了图像重建；实验证明，混合变差算法与本发明研制的装置相结合，可以得出分辨率较高的重建图像。

具体实施过程如下：

在计算机与成像算法模块1输入确定用于激励人体胸腔的安全电流的幅值和频率，通过电极阵列8将激励电流信号施加至人体胸腔表面，对人体胸腔即被测目标10进行激励，通过电极阵列8测量胸腔表面感应的电压信号并送入计算机，在计算机内部通过混合变差算法进行成像，总控模块3为系统控制核心，接受来自计算机和成像算法模块1的控制指令，实现对生物电阻抗成像技术硬件系统的整体控制和协调工作，信号发生模块4产生频率范围是1kHz-1MHz范围内可调的正弦波信号；总控模块3控制整个系统的功能：在电荷泵电流源模块5将电压信号转换成电流信号，电流信号的幅值为0.1mA-1mA可调，通过激励通道选通模块6实现激励电极的选通与关闭，使电流信号按照设定的方式，流入相应的电极，注入被测目标10，通过控制测量通道选通模块7选通电极阵列8中的相应电极，提取被测目标表面感应出的电压信号，以便将信号送入信号调理模块9，控制信号调理模块9内部的可编程增益放大电路的放大倍数，使得信号调理模块9将感应电压信号进行有效放大后，再由信号调理模块9内部的滤波电路和乘法解调电路分别对感应的电压信号进行滤波和解调，并由信号调理模块9内部的模数转换电路将其转化为数字量信号；然后再将转化的数字量电压信号由通信模块2送入计算机与成像算法模块1，将数字量电压信号复原为模拟量电压信号，通过混合变差算法实现胸腔电阻抗图像重建，最后由计算机输出重建胸腔电阻抗图像。

上述实施例中进行混合变差生物电阻抗成像所用装置是通过公知途径获得的，所涉及的操作方法是本领域技术人员所能掌握的。

Claims (2)

1.一种混合变差生物电阻抗成像方法，其特征在于步骤如下：

A.设置进行混合变差生物电阻抗成像所用装置：

进行混合变差生物电阻抗成像所用装置的构成包括计算机与成像算法模块、通信模块、总控模块、信号发生模块、电荷泵电流源模块、激励通道选通模块、测量通道选通模块、电极阵列和信号调理模块；计算机与成像算法模块用来控制总控模块和程序运行，总控模块通过控制总线、数据总线和地址总线完成对生物电阻抗成像技术系统各部分精确实时控制和快速信息的交换，信号发生模块用于产生正弦电压信号，电荷泵电流源模块将电压转换成电流，激励通道选通模块用于选择和切换电极阵列中的不同电极注入模式以将电流信号注入被测目标，测量通道选通模块用于选择和切换电极阵列的测量模式，以提取被测目标表面感应的电压信号，送至信号调理模块，信号调理模块通过其内部的可编程增益放大电路、滤波电路、乘法解调电路和模数转换电路，对电压信号进行放大、滤波、解调和模数转换，转换为数字量的电压信号由通信模块送入总控模块，并由通信模块实时地送入计算机与成像算法模块，在计算机内调用成像算法完成图像的重建工作。

B.进行混合变差生物电阻抗成像方法以实现对胸腔的电阻抗图像重建，其技术方案如下：

第一步，正问题求解：

(1.1)建立基于有限元方法的生物电阻抗成像技术的数值模型(以下简称数值模型)；

(1.2)定义数值模型的物理特性；

(1.3)对数值模型进行离散；

(1.4)施加边界条件；

(1.5)计算边界电压；

第二步，确定算法的目标函数：

(2.1)确定吉洪诺夫算法的目标函数；

(2.2)确定总变差算法的目标函数；

(2.3)确定混合变差算法的目标函数；

第三步，进行逆问题求解和胸腔电阻抗图像重建：

首先进行逆问题求解，即采用最速下降法求解混合变差算法的目标函数；

然后，包括采用以下的胸腔电阻抗图像重建算法，

(3.1)采用吉洪诺夫算法进行胸腔电阻抗图像重建，

(3.2)采用总变差算法进行胸腔电阻抗图像重建，

(3.3)采用混合变差算法进行胸腔电阻抗图像重建，步骤如下：

(3.3.1)临床测得电压信号或实验仿体模型测得电压信号或由仿真计算得到电压数据；

(3.3.2)根据混合变差目标函数，应用优化的L曲线方法，自适应调节正则化参数λ²；

(3.3.3)重建胸腔电导率分布，获得重建的胸腔电阻抗图像；

第四步，输出图像：

由上述A所述的计算机输出重建胸腔电阻抗图像；

C.用上述A所述的进行混合变差生物电阻抗成像所用装置完成上述B所述的进行混合变差生物电阻抗成像方法以实现对胸腔的图像重建的技术方案的实施过程如下：

在计算机与成像算法模块输入确定用于激励人体胸腔的安全电流的幅值和频率，通过电极阵列将激励电流信号施加至人体胸腔表面对人体胸腔即被测目标进行激励，通过电极阵列测量胸腔表面感应的电压信号并送入计算机，在计算机内部，通过混合变差算法进行成像，总控模块为系统控制核心，接受来自计算机和成像算法模块的控制指令，实现对生物电阻抗成像技术硬件系统的整体控制和协调工作，信号发生模块产生1kHz-1MHz范围内可调的正弦波信号；总控模块控制整个系统的功能：在电荷泵电流源模块将电压信号转换成电流信号，电流信号的幅值位于0.1mA至1mA之间且可调节，通过激励通道选通模块，实现激励电极的选通与关闭，使电流信号按照设定的方式流入相应的电极，注入被测目标；通过控制测量通道选通模块选通电极阵列中的相应电极提取被测目标表面感应出的电压信号，以便将信号送入信号调理模块；控制信号调理模块内部的可编程增益放大电路的放大倍数，使得信号调理模块将感应电压信号进行有效放大后，再由信号调理模块内部的滤波电路和乘法解调电路分别对感应的电压信号进行滤波和解调，并由信号调理模块内部的模数转换电路将其转化为数字量信号；然后再将转化的数字量电压信号由通信模块送入计算机与成像算法模块，将数字量电压信号复原为模拟量电压信号，通过混合变差算法实现胸腔电阻抗图像重建，最后由计算机输出重建胸腔电阻抗图像。

2.根据权利要求1所述一种混合变差生物电阻抗成像方法，其特征在于：所述根据混合变差目标函数，应用优化的L曲线方法，自适应调节正则化参数λ²，其具体的方法如下：

在优化的L曲线中，横坐标ξ是残差，纵坐标η是正则解，即

(ξ(λ²),η(λ²))＝(log||y-Ax||₂,log||x||₂)，

其中y为电压分布，A为雅克比矩阵，x为目标电导率分布；

对于混合变差正则化，优化的L曲线具有三个特征不同的部分，上部、中部和下部，其中上部和下部几乎垂直，中部几乎水平，优化的L曲线有两个拐角，其中，上部“垂直”部分对应着正则化参数过小时的解，下部“垂直”部分对应着正则化参数过大时的解，中部“水平”部分则对应着正则化参数取值适中时的解；优化的L曲线具有两个拐角，两个拐角将曲线的“水平”部分与上部“垂直”部分、下部“垂直”部分分开，通过定位两个拐角即可确定优化的L曲线的“水平”部分的范围，此范围对应着正则化误差和拟合误差之间达到平衡时的解；

上述定位拐角的方法是依据曲率最值准则，即选择优化的L曲线上曲率最大值点和最小值点，这两个拐角由下式求出：

$\mathrm{\&kappa;}\left({\mathrm{\&lambda;}}^2\right)=\arg max\frac{{\mathrm{\&xi;}}^{\&prime;}{\mathrm{\&eta;}}^{\&prime;\&prime;}-{\mathrm{\&xi;}}^{\&prime;\&prime;}{\mathrm{\&eta;}}^{\&prime;}}{{({\left({\mathrm{\&xi;}}^{\&prime;}\right)}^2+{\left({\mathrm{\&eta;}}^{\&prime;}\right)}^2)}^{3/2}}$

其中ξ'、η'表示ξ、η的一阶导数，ξ”、η”表示ξ、η的二阶导数，

这样，由优化的L曲线“水平”部分来确定合适的正则化参数，

具体步骤如下：

(1)赋初值λ²＝10^-6，a＝0，计算混合变差算法目标函数的最小值，得到逆问题求解结果；

(2)将第一步的结果作为初值，进行下一步的迭代，根据第一步重构的结果，分别将a依次赋值为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1，针对每个给定的a值，依次寻找其对应的正则化参数λ²，并计算此时混合变差算法目标函数的最小值，重建电导率分布；当满足下式的停止准则时，迭代停止；

$\sqrt{\frac{\underset{i=1}{\overset{m}{\&Sigma;}}{\&lsqb;x\left(i\right)-x_0\left(i\right)\&rsqb;}^2}{\underset{i=1}{\overset{m}{\&Sigma;}}{\&lsqb;x\left(i\right)\&rsqb;}^2}}<{\mathrm{\&tau;\&delta;}}_n,$

其中，x(i)表示重建的电导率分布，x₀(i)表示设定的初始电导率分布，m表示电导率分布矢量的维数，τ为1.0至1.5之间的一个预设常数，δ_n为电压数据的噪声水平；

在利用最速下降法对混合变差算法的目标函数进行求解过程中，由于总变差算法的目标函数不是光滑函数，并非处处可导，本发明中为了解决时不可微的情况，采用如下近似：

$\left|\right|R\&CenterDot;\&dtri;x|{|}_1\&ap;{\&Integral;}_{\mathrm{\&Omega;}}\sqrt{{(R\&CenterDot;\&dtri;x)}^2+\mathrm{\&beta;}},$

其中R是正则化矩阵，Ω为求解场域，是梯度算子，β为一个选定的常数，则有：

$\begin{array}{c}{F}^{\&prime;}\left(x\right)=\{{(Ax-y)}^T(Ax-y)+{\mathrm{\&lambda;}}^2(1-a){\&lsqb;R(x-x_0)\&rsqb;}^T\&lsqb;R(x-x_0)\&rsqb;+{\mathrm{\&lambda;}}^2\&CenterDot;a{\&Integral;}_{\mathrm{\&Omega;}}\sqrt{{(R\&CenterDot;\&dtri;x)}^2+\mathrm{\&beta;}}d\mathrm{\&Omega;}\}\\ =x^T{A}^{T}A-y^{T}A+{\mathrm{\&lambda;}}^2\&lsqb;(1-a){(x-x_0)}^T\&CenterDot;R^T\&CenterDot;R+a\&dtri;\&CenterDot;\frac{R\&CenterDot;\&dtri;x}{\sqrt{|\&dtri;x{|}^2}+\mathrm{\&beta;}}\&rsqb;\end{array},$

利用最速下降法求出混合变差算法的目标函数，求解中所依据的公式为：

x^k+1＝x^k-tF'(x)，

其中，t为最速下降法的步长，k是当前迭代的次数，每次迭代中采用线性搜索方法，使得不等式F^k+1(x)＜F^k(x)成立。