CN109662712A - 一种基于阻抗频谱约束的动态多频电阻抗断层成像方法 - Google Patents

Abstract

本发明属于电阻抗断层成像领域，具体涉及一种基于阻抗频谱约束的动态多频电阻抗断层成像方法,步骤1：在成像区域内进行有限面元剖分，确定成像域内的组织种类数和每种组织的频谱，建立体积分数和边界电压的关系模型；步骤2：采用不同频率的激励电流对成像区域边界处进行激励，测量每个频率下不同时刻的边界电压数据；步骤3：根据步骤1得到的体积分数变化和边界电压变化的关系构造逆问题目标函数并求解，得到成像区域内体积分数变化值；步骤4：将体积分数变化转化为阻抗变化并利用色阶表示在成像区域内，得到阻抗变化图像。本发明从增加观测数据量的角度提供了一种改善逆问题病态性的新方案，进而提高了重构图像的空间分辨率和抗噪能力。

Description

一种基于阻抗频谱约束的动态多频电阻抗断层成像方法

技术领域

本发明属于电阻抗断层成像领域，具体涉及一种基于阻抗频谱约束的动态 多频电阻抗断层成像方法。

背景技术

电阻抗断层成像技术通过安装在物体表面的电极对物体有规律地施加激 励电流，若物体内部存在阻抗变化，则会引起表面测量电极电压的变化。然后 基于有限元分析建立正问题模型，得到边界电压变化同成像域内阻抗变化的对 应关系，最后根据测量得到的边界电压变化结合相应重建算法可以获得物体内 部阻抗变化图像。在实际应用中，病态性严重是目前各种成像方法面临的主要 问题，病态性的实质是物体内部阻抗分布和表面测量电压之间存在非线性关系， 且为了获得数值解，需要对这种非线性关系进行离散化处理，而离散化处理存 在着一定的误差，这将导致实际测量电压即使存在一个很小的噪声都会导致内 部重构阻抗的剧烈变化，这非常不利于我们实现临床上连续动态的EIT监测。

申请号为：CN105232044A，名称为《一种自构造背景帧的电阻抗断层成像 方法》的专利提出利用动态成像对两帧测量数据求差来削弱噪声的影响以改善 病态性，但是该方法测量的数据远远小于待求数据，无法从根源上改善病态性， 图像的抗噪性能和空间分辨率较差。

申请号为：CN102894961A，名称为《一种多频电阻抗断层成像的谱成像方 法》的专利提出利用多频信息来做独立成分分析，但是该方法重构的是内部电 导率随频率的变化量Δδ并且是一个静态时刻下的频差成像，仍然无法从根源 上改善病态性，还不利于对疾病的实时监测。

发明内容

针对现有技术中存在的无法从根本上改善病态性问题和图像质量差的问 题，本发明提供了一种基于阻抗频谱约束的动态多频电阻抗断层成像方法，包 括如下步骤：

步骤1：在成像区域内进行有限面元剖分，确定成像域内的组织种类数， 获得每种组织的体积分数向量和阻抗频谱以及体积分数和阻抗的线性关系，建 立体积分数和边界电压的关系模型；

步骤2：采用不同频率的激励电流对成像区域边界处进行激励，测量每个 频率下不同时刻的边界电压数据，采用式(1)获得每个频率下的一组边界电 压差数据：

Δυ_i(t)＝υ_i(t)-υ_i(1) (式1)，

其中，υ_i(t)表示激励频率为ω_i的不同时刻下的边界电压数据，i表示阻抗 频谱内的第i个频率，i≤I且I为正整数，t≥2，υ_i(1)表示激励频率为ω_i的时 刻1下边界电压数据；

步骤3：根据步骤1的体积分数和边界电压的关系模型构造逆问题目标函数， 优化目标函数并利用逆问题方法得到体积分数变化的求解公式，将步骤2得到 的多个频率的边界电压差数据输入求解公式，得到成像区域内体积分数变化值；

步骤4：将步骤3得到的成像区域内体积分数变化值代入步骤1中得到的体 积分数和阻抗的线性关系，获得成像区域内的阻抗变化，将阻抗变化利用色阶 表示在成像区域内，得到阻抗变化图像。

进一步的，步骤1包括如下子步骤：

步骤1.1：在成像区域内进行有限面元剖分，得到多个面元，确定成像域 内的组织种类数，并设定每个面元内每种组织的体积分数值，获得体积分数矩 阵，将体积分数矩阵向量化得到体积分数模型；

步骤1.2：获得每种组织的阻抗频谱，结合步骤1.1获得的体积分数模型， 获得体积分数和阻抗的线性关系，如式2所示：

σ(ω_i)＝A_iF (式2)

其中，σ是阻抗，ω_i为阻抗频谱中第i个激励频率，A_i是ω_i下的阻抗频谱构 成的系数矩阵，F是体积分数向量；

步骤1.3：将步骤1.2得到的体积分数和阻抗的关系结合阻抗变化和边界电 压变化的关系，得到如式3的体积分数变化和边界电压变化的关系：

J_iA_iΔF＝Δυ_i (式3)

其中，J_i是ω_i下的雅克比矩阵，ΔF是体积分数向量的变化量，Δυ_i是ω_i下 的边界电压的变化量

一种基于阻抗频谱约束的动态多频电阻抗断层成像系统，包括多频数据采 集模块、体积分数构造模块、图像重构模块和图像输出模块；

所述多频数据采集模块包括电极、导联线、多路开关、频率控制器和通信 处理器；

所述体积分数构造模块包括体积分数模型器和通信处理器，该模块基于有 限元剖分构建待测体的体积分数模型，并向计算机输出该模型下体积分数变化 和边界电压变化之间的关系；

所述图像重构模块包括目标函数构造优化器、目标函数求解器、边界条件 限制器以及通信处理器，该模块利用体积分数构造模块输出的体积分数变化和 边界电压变化之间的关系构造一个目标函数并将之优化，然后将从计算机接收 的边界电压数据输入目标函数求解器求解体积分数变化，最后经过边界条件限 制，得到最终的体积分数变化值；

所述图像输出模块包括体积分数阻抗转化器、画图器、显示器以及通信处 理器，该模块用于接收计算机输出的体积分数变化并利用转化器将其转化成阻 抗变化并画图显示。

本发明具有以下技术效果：

(2)本发明从增加观测数据量的角度提供了一种改善逆问题病态性的新 方案，进而提高了重构图像的空间分辨率和抗噪能力；首先该方法根据组织的 频谱特性，构建体积分数模型，其次建立体积分数和边界电压之间的关系，最 终同时利用多个激励频率下的测量电压差数据重构一帧时差图像，以提高重建 图像空间分辨率并减少噪声带来的图像伪影。

(2)本发明的方案经仿真测试验证，对求解域内目标位置定位较传统阻 抗成像更准确，且具有更强的噪声抑制能力，有望对多频动态电阻抗断层成像 提供一个新的发展方向。

附图说明

图1是本发明的系统连接图；

图2是圆域有限元剖分模型和电极位置示意图；

图3是距圆域中心不同距离的四个目标示意图；

图4成像结果对比图；

图5是两组织的阻抗频谱特性；

图6是阻抗频谱构成的系数矩阵。

具体实施方式

一种基于阻抗频谱约束的动态多频电阻抗断层成像方法，包括如下步骤：

步骤1：在成像区域内进行有限面元剖分，确定成像域内的组织种类数， 获得每种组织的体积分数向量和阻抗频谱以及体积分数和阻抗的线性关系，建 立体积分数和边界电压的关系模型；成像区域一般指待观测的人体，具体取决 于成像目标所在部位，如果是肺部EIT，成像区域即肺部一横截面，具体即电 极放置的横截面，如果是脑部EIT，成像区域即脑部一横截面，具体即电极放 置的横截面；

步骤2：采用不同频率的激励电流对成像区域边界处进行激励，测量每个 频率下不同时刻的边界电压数据，采用式(1)获得每个频率下的一组边界电 压差数据：

Δυ_i(t)＝υ_i(t)-υ_i(1) (式1)，

其中，υ_i(t)表示激励频率为ω_i的不同时刻下的边界电压数据，i表示阻抗 频谱内的第i个频率，i≤I且I为正整数，t≥2，υ_i(1)表示激励频率为ω_i的时 刻1下边界电压数据；

其中，成像区域边界一般指待观测的人体表面，例如肺部皮肤或者脑部头 皮；

步骤3：根据步骤1得到的体积分数变化和边界电压变化的关系构造逆问题 目标函数，使得重构误差的二范数最小，优化目标函数并利用逆问题方法得到 体积分数变化的求解公式，将步骤2得到的多组边界电压差输入求解公式，得 到成像区域内体积分数变化值；

步骤4：将步骤3得到的成像区域内体积分数变化值代入步骤1中得到的体 积分数和阻抗的线性关系，获得成像区域内的阻抗变化，将阻抗变化利用色阶 表示在成像区域内，得到阻抗变化图像。

具体的，步骤1包括如下子步骤：

步骤1.1：在成像区域内进行有限面元剖分，得到多个面元，确定成像域 内的组织种类数，并设定每个面元内每种组织的体积分数值，获得体积分数矩 阵，将体积分数矩阵向量化得到体积分数模型；

步骤1.2：，获得每种组织的阻抗频谱，由于混合组分的阻抗可以表示成各 组分阻抗加权和的形式，结合步骤1.1获得的体积分数模型，获得体积分数和 阻抗的线性关系，如式2所示：

σ(ω_i)＝A_iF (式2)

其中，σ是阻抗，ω_i为阻抗频谱中第i个激励频率，A_i是ω_i下的阻抗频谱 构成的系数矩阵如图6所示

F是体积分数向量，且F∈R^(T*N)*1，其形式如下：

其中，f_nj为面元n内第j种组织的体积分数值且n∈[1，N]，j∈[1，T]；

步骤1.3：将步骤1.2得到的体积分数和阻抗的关系结合阻抗变化和边界电 压变化的关系，得到如式3的体积分数变化和边界电压变化的关系，该关系是 不含噪声理想状态下的体积分数模型的离散正问题：

J_iA_iΔF＝Δυ_i (式3)

其中，J_i是ω_i下的雅克比矩阵，ΔF是体积分数向量的变化量，Δυ_i是ω_i下 的边界电压的变化量。

优选的，成像种类包括：脑部电阻抗成像、肺部电阻抗成像，脑部电阻抗 成像包括的组织有：正常脑实质，缺血脑组织，出血脑组织，肺部电阻抗成像 包括的组织有：正常肺组织，充气肺组织，放气肺组织。

优选的，我们已知在传统阻抗重构时，不含噪声理想状态下的传统的离散 正问题即为所述阻抗和边界电压的关系：

J_iΔσ_i＝Δv_i。

具体的，步骤2在成像区域边界水平放置16个电极，选择临近激励和对侧 测量模式，施加多个不同频率的激励频率，可得到多个频率下的边界电压数据， 再选取某一个时刻的边界电压数据作为背景帧，将任意时刻的边界电压数据和 背景帧作差可得到该时刻的边界电压数据的变化值。

具体的，步骤3包括如下子步骤：

步骤3.1：根据步骤1得到的体积分数变化和边界电压变化的关系构造目标 函数Φ使得重构误差的二范数最小，为了将I组边界电压差同时应用于一帧时 差成像，使得所有频率下的数据差二范数最小，目标函数Φ为：

其中，λ是正则化参数，R是正则化矩阵,ΔF体积分数变化向量，Δυ是任意 时刻和背景帧之间的边界电压差向量，S是组装矩阵内含有I个元素矩阵，每 个元素矩阵由每个频率下的雅克比矩阵和系数矩阵相乘得到，S形如：

步骤3.2：对目标函数Φ进行优化，使其满足：任何一个面元内所有体积 分数和为1，即并且任何一个体积分数值在[0,1]之间，即 得到优化后的目标函数Φ′；

步骤3.3：利用逆问题求解方法求解优化后的目标函数Φ′得到求解式，将

步骤2采集到的边界电压差Δυ输入求解式中，得到成像区域内T-1种组织的体 积分数变化ΔF_T-1。

优选的，步骤3.3具体包括如下子步骤：

步骤3.3.1：设定体积分数初始分布

步骤3.3.2：采取Gauss-Newton方法得到目标函数的步进方向d；

步骤3.3.3：采取Brent线性搜索法寻找最佳的步进长度β；

步骤3.3.4：利用步骤3.3.2和步骤3.3.3得到的d和步进长度β求得当前时 刻计算体积分数分布令中每个元素大于等于0 小于等于1，得到当前时刻真实体积分数分布F_T-1；

步骤3.3.5：利用F_T-1减去得到区域内部体积分数变化量ΔF_T-1。

一种基于阻抗频谱约束的动态多频电阻抗断层成像系统，包括多频数据采 集模块、体积分数构造模块、图像重构模块和图像输出模块。

所述多频数据采集模块包括电极、导联线、多路开关、频率控制器和通信 处理器；该模块按照计算机的指令利用导线向激励电极施加激励电流，并利用 测量电极采集边界电压数据，多路开关负责切换激励电极和测量电极而频率控 制器负责切换激励频率，最后将采集的多频边界电压数据传输给计算机控制模 块。

所述体积分数构造模块包括体积分数模型器和通信处理器，该模块基于有 限元剖分构建待测体的体积分数模型并向计算机输出该模型下体积分数变化 和边界电压变化之间的关系。

所述图像重构模块包括目标函数构造优化器、目标函数求解器、边界条件 限制器以及通信处理器，该模块利用体积分数构造模块输出的体积分数变化和 边界电压变化之间的关系构造一个目标函数并将之优化，之后利将从计算机接 收的边界电压数据输入目标函数求解器求解体积分数变化，最后经过边界条件 限制，得到最终的体积分数变化值。

所述图像输出模块包括体积分数阻抗转化器、画图器、显示器以及通信处 理器，该模块接收计算机输出的体积分数变化并利用转化器将其转化成阻抗变 化，最后画图显示。

以下给出本发明的具体实施方式，需要说明的是本发明并不局限于以下具 体实施例，凡在本申请技术方案基础上做的等同变换均落入本发明的保护范围。

仿真方案：

将采用包含两组织的单目标仿真，按照如下步骤实施：

(1)如图1设置两种阻抗频谱不同的组织，并选取其中的三个频点来进行 仿真，其中t₁组织是背景，t₂组织是目标。

(2)建立圆域仿真模型，基于VS2013平台对圆域进行有限元剖分，并在 圆域边界均匀设置16个电极如图1，得到体积分数成像仿真模型。

(3)所有面元的f₁＝1；f₂＝0(即所有面元被t₁组织占据)，得到初始 化向量F⁰。装配各个频率下的系数矩阵A_i并分别乘以向量F⁰得到该频率下的 阻抗向量σ_i。由正问题映射得到各个频率下的边界电压数据υ_i，组成背景帧数 据。

(4)在背景基础上设置成像目标其f₁＝0；f₂＝1(即所有面元被t₂组织 占据)，同(3)的方法生成前景帧数据。为了较好的测试该成像方法性能，施 例选取了4种测试目标，它们距离圆域中心距离不同，如图3。将前景帧数据和 背景帧数据做差分得到边界电压的变化量。

(5)分别用传统的阻抗成像方法和基于频谱约束的多频阻抗成像方法对 上述4种目标成像。

(6)对比两种成像方法下的重构图像质量，具体见图4，其中第一列是目 标真值图，第二列是传统阻抗成像方式下的重构图，第三列是基于频谱约束的 多频阻抗成像方法下的重构图。

实施例：

应用在脑部EIT时，步骤1：如图2-图5所示，首先假设三维区域Ω内含有T 种组织，分别为：t₁……t_j……t_T，(j和T是正整数且j≤T)并获得这T种组织 的阻抗频谱，对Ω做有限元剖分，由于面元内的阻抗是连续一致的，得到阻抗 向量σ_i∈R^N*1，其中N是剖分面元数，i是激励频率，R^N*1是大小为N*1的矩阵, 在每个面元内设定每种组织的体积分数值，获得体积分数矩阵FM∈R^T*N，其 中，R^T*N为大小为T*N的矩阵,R^T*N的列是面元n的各组织占比,且n≤N，行 是组织j在各个面元的占比，之后将FM按列向量化得到体积分数向量 F＝vec(FM),F∈R^(T*N)*1；

步骤2：在该三维区域中心水平面边界上等间距放置16个电极，任意时刻 采用I个激励频率激励电流，分别表示为ω₁……ω_i……ω_I，其中，i、I均为 正整数且i≤I每个激励频率下可以得到一个边界电压向量υ_i∈R^K*1，其中K是 一个激励频率下的所有边界电压值个数，将这I个边界电压向量纵向拼接可得 该时刻的边界电压向量υ(t)，获得每个频率下的一组边界电压差数据如式(1) 所示：

Δυ_i(t)＝υ_i(t)-υ_i(t₀) (式1)，

其中υ_i(t)表示激励频率为ω_i的不同时刻下的边界电压数据，且t≥时刻2， υ_i(t₀)表示激励频率为ω_i的时刻1下边界电压数据；

步骤3.1：根据步骤1得到的体积分数变化和边界电压变化的关系构造目标 函数Φ使得重构误差的二范数最小，为了将I组边界电压差同时应用于一帧时 差成像，使得所有频率下的数据差二范数最小，目标函数Φ为：

其中，λ是正则化参数，R是正则化矩阵,ΔF体积分数变化向量，Δυ是任意 时刻和背景帧之间的边界电压差向量，S是组装矩阵内含有I个元素矩阵，每 个元素矩阵由每个频率下的雅克比矩阵和系数矩阵相乘得到，S形如：

步骤3.2：对目标函数Φ进行优化，在目标函数中替换使 其满足：任何一个面元内所有体积分数和为1，即对应的t₁组 织通常是背景组织，在替换f₁之后，未知数个数变小，所以待求的未知数由ΔF 变成ΔF_T-1，此时将重构T-1种组织(f₂……f_j……f_T)的体积分数变化量ΔF_T-1， 相应的系数矩阵也变为A′；同时要对迭代结果施加边界条件，以保证F_T-1中的 任何一个体积分数值在[0,1]之间，即得到优化后的目标函数Φ′；

步骤3.3：利用逆问题求解方法求解优化后的目标函数Φ′得到求解式，将 步骤2采集到的边界电压差Δυ输入求解式中，得到成像区域内T-1种组织的体 积分数变化ΔF_T-1。

步骤4：将步骤3得到的成像区域内体积分数变化值代入步骤1中得到的体 积分数和阻抗的线性关系，获得成像区域内的阻抗变化，将阻抗变化利用色阶 表示在成像区域内，得到阻抗变化图像，如图4所示。

Claims (3)

1.一种基于阻抗频谱约束的动态多频电阻抗断层成像方法，其特征在于，包括如下步骤：

步骤1：在成像区域内进行有限面元剖分，确定成像域内的组织种类数，获得每种组织的体积分数向量和阻抗频谱以及体积分数和阻抗的线性关系，建立体积分数和边界电压的关系模型；

步骤2：采用不同频率的激励电流对成像区域边界处进行激励，测量每个频率下不同时刻的边界电压数据，采用式(1)获得每个频率下的一组边界电压差数据：

Δυ_i(t)＝v_i(t)-v_i(1) (式1)，

其中，υ_i(t)表示激励频率为ω_i的不同时刻下的边界电压数据，i表示阻抗频谱内的第i个频率，i≤I且I为正整数，t≥2，υ_i(1)表示激励频率为ω_i的时刻1下边界电压数据；

步骤3：根据步骤1的体积分数和边界电压的关系模型构造逆问题目标函数，优化目标函数并利用逆问题方法得到体积分数变化的求解公式，将步骤2得到的多个频率的边界电压差数据输入求解公式，得到成像区域内体积分数变化值；

步骤4：将步骤3得到的成像区域内体积分数变化值代入步骤1中得到的体积分数和阻抗的线性关系，获得成像区域内的阻抗变化，将阻抗变化利用色阶表示在成像区域内，得到阻抗变化图像。

2.如权利要求1所述的基于阻抗频谱约束的动态多频电阻抗断层成像方法，其特征在于，步骤1包括如下子步骤：

步骤1.1：在成像区域内进行有限面元剖分，得到多个面元，确定成像域内的组织种类数，并设定每个面元内每种组织的体积分数值，获得体积分数矩阵，将体积分数矩阵向量化得到体积分数模型；

步骤1.2：获得每种组织的阻抗频谱，结合步骤1.1获得的体积分数模型，获得体积分数和阻抗的线性关系，如式2所示：

σ(ω_i)＝A_iF (式2)

其中，σ是阻抗，ω_i为阻抗频谱中第i个激励频率，A_i是ω_i下的阻抗频谱构成的系数矩阵，F是体积分数向量；

步骤1.3：将步骤1.2得到的体积分数和阻抗的关系结合阻抗变化和边界电压变化的关系，得到如式3的体积分数变化和边界电压变化的关系：

J_iA_iΔF＝Δυ_i (式3)

其中，J_i是ω_i下的雅克比矩阵，ΔF是体积分数向量的变化量，Δυ_i是ω_i下的边界电压的变化量。

3.一种基于阻抗频谱约束的动态多频电阻抗断层成像系统，其特征在于：包括多频数据采集模块、体积分数构造模块、图像重构模块和图像输出模块；

所述多频数据采集模块包括电极、导联线、多路开关、频率控制器和通信处理器；

所述体积分数构造模块包括体积分数模型器和通信处理器，该模块基于有限元剖分构建待测体的体积分数模型，并向计算机输出该模型下体积分数变化和边界电压变化之间的关系；

所述图像重构模块包括目标函数构造优化器、目标函数求解器、边界条件限制器以及通信处理器，该模块利用体积分数构造模块输出的体积分数变化和边界电压变化之间的关系构造一个目标函数并将之优化，然后将从计算机接收的边界电压数据输入目标函数求解器求解体积分数变化，最后经过边界条件限制，得到最终的体积分数变化值；

所述图像输出模块包括体积分数阻抗转化器、画图器、显示器以及通信处理器，该模块用于接收计算机输出的体积分数变化并利用转化器将其转化成阻抗变化并画图显示。