CN103065286A - 一种准静态电阻抗成像中的图像重构方法 - Google Patents

Abstract

本发明公开了一种准静态电阻抗成像中的图像重构方法，采用电阻抗成像方法，对待测目标分别在频率为和

的电流激励下进行激励，采集两次电流激励下的边界电压

和

Description

一种准静态电阻抗成像中的图像重构方法

技术领域

本发明属于电阻抗断层成像技术领域，涉及一种准静态电阻抗成像中的图像重构方法。

背景技术

电阻抗断层成像技术（electrical impedance tomography,EIT）是一种新型的医学成像技术，通过向人体施加安全的激励电流（电压），同时测量体表电压（电流）信号，经重构方法获得人体某一断层内部的阻抗或阻抗变化分布。

电阻抗断层成像可分为：静态成像（以阻抗的绝对值分布为成像目标）和动态成像（以阻抗的变化为成像目标）。由于静态成像对数据采集系统要求苛刻，在实际应用中受到了限制，目前仍然处于实验室研究阶段。动态成像通过对不同时刻的测量数据进行差分成像，一定程度上减少了接触阻抗等系统误差对图像重构的干扰，所以可用于对某些疾病的发展过程进行实时图像监测。但是对于病灶已形成、在短时间内不会发生阻抗动态变化的疾病而言，采用动态EIT成像的方法不能快速检测，因此有必要发展一种新的EIT成像方法。

准静态电阻抗成像技术是一种新颖的电阻抗成像模式，通过在同一时刻向被测体施加两种或多种频率的激励电流，利用不同频率下的测量数据进行差分成像，重建出反应被测体内部阻抗分布，达到对疾病检测的目的。

人体各组织都具有一定的阻抗频谱特性，所以不同激励频率下的测量数据中不仅包含了病变信息而且具有正常组织的信息，因而不同频率下的数据的差分结果中包含了大量的非病变组织阻抗的信息，因此，会导致重构图像中存在大量的伪影，甚至无法分辨出病变。但是，在临床应用中，医务工作者通常更加关注图像中病变的信息(例如，位置和大小)而不是正常组织。

因此，为了消除非病变组织对成像结果的影响，凸显病变组织，以满足临床实际需求，需要一种更加合理、有效的准静态电阻抗成像方法。

发明内容

本发明解决的问题在于提供一种准静态电阻抗成像中的图像重构方法，该方法能够在成像结果中凸显病变组织（扰动目标），减少或消除正常组织（非目标区域）的影响。

本发明是通过以下技术方案来实现：

一种准静态电阻抗成像中的图像重构方法，包括以下步骤：

1）采用电阻抗成像方法，对待测目标分别在频率为

和

的电流激励下进行激励，采集两次电流激励下的边界电压

和

然后对所采集的数据进行加权求差处理，得到不同频率激励下的目标区域阻抗变化差异所引起的边界电压变化ΔV：

$\mathrm{\ΔV}=a\·V_{{\mathrm{\ω}}_2}^{\′}-V_{{\mathrm{\ω}}_1}^{\′}$

其中，a为两种不同频率下成像区域内阻抗的变化比例系数，

n为电极数；m为一次激励所产生的边界测量电压数；

2）然后计算不同频率激励下的目标区域阻抗分布变化Δρ：

其中，

为激励频率为

对应的敏感性矩阵；

利用不同频率激励下的目标区域阻抗分布变化Δρ重建出EIT图像。

还根据检测对象估计其场域内阻抗分布为ρ_α，则Δρ为：

$\mathrm{\Δ\ρ}=a\·{\left(\frac{{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}}{{\mathrm{\ρ}}_{\mathrm{\α}}}\right)}^2{S}^{-1}\·\mathrm{\ΔV}$

其中，

和

相等或接近，real(ρ)和imag(ρ)表示ρ的实部和虚部，

为实常数，敏感矩阵S的元素通过场域内阻抗分布ρ_α计算：

所述的ρ_α的值与实际检测对象的阻抗分布值成比例，使

为实常数，或是虚部与实部相比较很小，不影响成像结果。

所述在利用Δρ重建出EIT图像时，还需采用正则化方法加以修正。

所述采用的Tikhonov正则化方法加以修正：

$\mathrm{\Δ\ρ}=a\·{\left(\frac{{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}}{{\mathrm{\ρ}}_{\mathrm{\α}}}\right)}^2\·{(S^{T}S+\mathrm{\λR})}^{-1}{S}^T\·\mathrm{\ΔV},$ 其中，R为正则化矩阵，λ为正则化系数。

所述的

低于

为低频（通常选择

小于10KHz），

为高频（通常选择

大于10KHz）。

与现有技术相比，本发明具有以下有益的技术效果：

本发明提供的准静态电阻抗成像中的图像重构方法，在利用两组不同频率的测量数据成像前，考虑到成像区域内非目标信息对重构结果的影响，提出采用加权求差的数据处理方式，可以有效消除非目标信息在重构图像中产生的伪影，提高目标辨识度，提高重建图像质量。

为了与传统重建算法相比较，本发明提供的准静态电阻抗成像中的图像重构方法的重建结果清晰，而传统的图像重建算法的重建结果中存在大量伪影，甚至淹没了目标。

本发明提供的准静态电阻抗成像中的图像重构方法，具有一定的抗噪声能力，在仿真实验中添加的高斯白噪声不影响成像结果的清晰性。

附图说明

图1为有限元模型的示意图；

图2为经加权求差预处理后的电阻抗成像结果；

图3为传统的图像重建算法的电阻抗成像结果。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明提供的准静态电阻抗成像中的图像重构方法，包括以下步骤：

1）采用电阻抗成像方法，对待测目标分别在频率为

和的电流激励下进行激励，采集两次电流激励下的边界电压

和

然后对所采集的数据进行加权求差处理，得到不同频率激励下的目标区域阻抗变化差异所引起的边界电压变化ΔV；

2）然后计算不同频率激励下的目标区域阻抗分布变化Δρ，利用不同频率激励下的目标区域阻抗分布变化Δρ重建出EIT图像。

进一步，说明如下：

EIT的数学模型可由椭圆偏微分方程和混和边界条件给出：

是场域Ω内部的电势分布；I_l为注入电流；n为外法向单位矢量；ρ为场域内阻抗的分布。

EIT边界电压V和阻抗分布ρ的线性敏感性关系可表示为：

V=S·ρ       (2)

其中S为敏感矩阵，其元素为：

其中，i表示第i次的电流激励和电压测量，j表示第j个剖分单元，它是对一个单元的体积分；带有w的量表示区域内为均匀电导率分布在频率w下的量，为通过第i对激励测量电极对中激励电极输入电流后，在各单元中产生的电位梯度；为通过第i对激励测量电极对中测量电极输入电流后，在各单元中产生的梯度。

假设成像区域内阻抗均匀分布，在两种频率为

和

的电流激励下的阻抗分布分别为

和

相应的边界测量电压为

和

所以：

$V_{{\mathrm{\ω}}_1}=S_{{\mathrm{\ω}}_1}\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}$ $V_{{\mathrm{\ω}}_2}=S_{{\mathrm{\ω}}_2}\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_2}---\left(4\right)$

设 ${\mathrm{\ρ}}_{{\mathrm{\ω}}_2}=a\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_1},$ 根据（1）可知： $V_{{\mathrm{\ω}}_1}=a\·V_{{\mathrm{\ω}}_2}---\left(5\right)$

其中，a为两种不同频率下成像区域内阻抗的变化比例系数，

为频率为

的敏感性矩阵，

为频率为的敏感性矩阵。

如果成像区域内有一扰动目标，其阻抗频谱特性不同于背景，在频率

和的电流激励下，成像区域的整体阻抗分布为

和

边界电压为

与

由于目标的阻抗相较于整体场域内非目标的阻抗对敏感矩阵的影响很小，所以边界电压和阻抗分布关系可近似表示为：

$V_{{\mathrm{\ω}}_1}^{\′}=S_{{\mathrm{\ω}}_1}\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}^{\′}$ $V_{{\mathrm{\ω}}_2}^{\′}=S_{{\mathrm{\ω}}_2}\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_2}^{\′}---\left(6\right)$

设 $\mathrm{\Δ\ρ}={\mathrm{\ρ}}_{{\mathrm{\ω}}_2}^{\′}-a\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}^{\′}$ $\mathrm{\ΔV}=a\·V_{{\mathrm{\ω}}_2}^{\′}-V_{{\mathrm{\ω}}_1}^{\′}---\left(7\right)$

结合（4）（6）可知：非目标区域内Δρ≈0，目标区域内Δρ≠0；所以ΔV仅为不同频率激励下的目标区域阻抗变化差异所引起。

$\mathrm{\Δ\ρ}={\mathrm{\ρ}}_{{\mathrm{\ω}}_2}^{\′}-a\·{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}^{\′}=S_{{\mathrm{\ω}}_2}^{-1}\·V_{{\mathrm{\ω}}_2}^{\′}-a\·S_{{\mathrm{\ω}}_1}^{-1}\·V_{{\mathrm{\ω}}_1}^{\′}---\left(8\right)$

由（2）（3）可知 $S_{{\mathrm{\ω}}_1}=a^2\·S_{{\mathrm{\ω}}_2}$ $S_{{\mathrm{\ω}}_2}^{-1}=a^2\·S_{{\mathrm{\ω}}_1}^{-1}$

所以： $\mathrm{\Δ\ρ}=a^2\·S_{{\mathrm{\ω}}_1}^{-1}\·V_{{\mathrm{\ω}}_2}^{\′}-a\·S_{{\mathrm{\ω}}_1}^{-1}\·V_{{\mathrm{\ω}}_1}^{\′}=a\·S_{{\mathrm{\ω}}_1}^{-1}\·(a\·V_{{\mathrm{\ω}}_2}^{\′}-\·V_{{\mathrm{\ω}}_1}^{\′})=a\·S_{{\mathrm{\ω}}_1}^{-1}\·\mathrm{\ΔV}$

即： $\mathrm{\Δ\ρ}=a\·S_{{\mathrm{\ω}}_1}^{-1}\·\mathrm{\ΔV}---\left(9\right)$

进一步，由(9)式可见，欲求Δρ需具有激励频率为

时，场域内阻抗分布

的先验信息，以计算

但是在实际中很难获得准确的

所以,为了克服此问题，假设场域内阻抗分布ρ_α，选择ρ_α时，尽量使

和

相等或接近，其中，real(ρ)和imag(ρ)表示ρ的实部和虚部，在实际应用中，可以提前根据检测对象估计出ρ_α，将(9)式进行改写为如下形式：

$\mathrm{\Δ\ρ}=a\·{\left(\frac{{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}}{{\mathrm{\ρ}}_{\mathrm{\α}}}\right)}^2{S}^{-1}\·\mathrm{\ΔV}---\left(10\right)$

其中，敏感矩阵的元素通过场域内阻抗分布ρ_α计算，即：因为

和相等或接近，故

为实常数，不影响成像结果；

由于目标阻抗变化相对整个成像区域较小，a通过整体电压边界确定。最小二乘估计是一种参数估计方法，该方法简单，通过最小化误差的平方和找到一组数据的最佳函数匹配。基于最小二乘原理，选择代价函数 $J=\underset{x=1}{\overset{n}{\mathrm{\Σ}}}\underset{y=1}{\overset{m}{\mathrm{\Σ}}}{(a\·V_{{\mathrm{\ω}}_2}^{\′}(x\·y)-V_{{\mathrm{\ω}}_1}^{\′}(x\·y))}^2,$ 通过求解，a可表示为：

$a=\frac{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}\underset{y=1}{\overset{m}{\mathrm{\Σ}}}(V_{{\mathrm{\ω}}_1}^{\′}(x\·y)\·V_{{\mathrm{\ω}}_2}^{\′}(x\·y))}{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}\underset{y=1}{\overset{m}{\mathrm{\Σ}}}{\left(V_{{\mathrm{\ω}}_2}^{\′}(x\·y)\right)}^2}---\left(11\right)$

其中，n为电极数；m为一次激励所产生的边界测量电压数；

然后采用合适的逆问题求解方法，通过(10)可计算出Δρ。

下面通过建立有限元模型来对上述方法的成像结果进一步说明。具体包括：

1、建立有限元模型，设置模型内阻抗分布。有限元模型如图2所示，在模型边界处等距离放置16个电极，模型半径为8cm，采用三角单元剖分，共512个剖分单元，采用对向激励和相邻测量的模式，激励电流为1mA。

当激励频率为f₁时，非目标区域的阻抗设置为0.5+0.05i（S/m），目标区域的电导率设置为0.15+0.01i（S/m）；

当激励频率为f₂时，非目标区域电阻设置为抗0.55+0.06i（S/m），目标区域的电导率设置为0.9+0.1i（S/m）。

以模拟待测对象在不同激励频率下阻抗分布的变化。

由于EIT对噪声非常敏感，因此为了更加接近实际，分别向仿真数据的实部和虚部加入0.1%的高斯白噪声。

2、计算敏感矩阵。在本实施例中，选择ρ_α=1+0.1i计算敏感矩阵。敏感矩阵的元素为：

i表示第i次的电流激励和电压测量，j表示第j个剖分单元，

为通过第i对激励测量电极对中激励电极输入电流后，在各单元中产生的电位梯度；

为通过第i对激励测量电极对中测量电极输入电流后，在各单元中产生的梯度。因为采用的激励测量模式为：对向激励，相邻测量，所以S矩阵大小为256×512。

$S=\left[\begin{array}{cccc}0& \·\·\·& \·\·\·& 0\\ -3.067e-04+6.196e-05i& \·\·\·& \·\·\·& -0.01280+0.002587i\\ -4.152e-04+8.389e-05i& \·\·\·& \·\·\·& -0.004157+0.0008398i\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ 4.023e-04-8.129e-05i& \·\·\·& \·\·\·& 5.107e-04-1.032e-04i\\ 3.203e-04-6.471e-05i& \·\·\·& \·\·\·& 9.836e-04-1.987e-04i\\ 0& \·\·\·& \·\·\·& 0\end{array}\right]$

3、模型设置完成后，对待测目标分别在频率为

和

的电流下进行激励，

为低频（通常选择

小于10KHz），为高频（通常选择

大于10KHz）；采用有限元方法，分别计算两种频率电流激励下的阻抗分布的边界电压V₁′和V₂′；

而在实际检测时中，V₁′和V₂′即为两种不同激励频率下的测量数据

和 $V_{{\mathrm{\ω}}_2}^{\′}.$

$V_{{\mathrm{\ω}}_1}^{\′}=\left[\begin{array}{c}0\\ -0.4707+0.04699i\\ -0.3046+0.03039i\\ -0.2565+0.02561i\\ -0.2558+0.02555i\\ \·\\ \·\\ \·\\ \·\\ \·\\ \·\\ 0.2556-0.02554i\\ 0.2561-0.02558i\\ 0.3048-0.03043i\\ 0.4690-0.04684i\\ 0\end{array}\right]$ $V_{{\mathrm{\ω}}_2}^{\′}=\left[\begin{array}{c}2.029e-04-3.020e-05i\\ -0.4226+0.04610i\\ -0.2730+0.02974i\\ -0.2295+0.02515i\\ -0.2299+0.02513i\\ \·\\ \·\\ \·\\ \·\\ \·\\ \·\\ 0.2304-0.02508i\\ 0.2298-0.02506i\\ 0.2736-0.02994i\\ 0.4219-0.04604i\\ 0.001034-7.990e-06i\end{array}\right]$

4、求解加权因子，并对数据采用加权求差的处理方式

利用 $a=\frac{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}\underset{y=1}{\overset{m}{\mathrm{\Σ}}}(V_{{\mathrm{\ω}}_1}^{\′}(x\·y)\·V_{{\mathrm{\ω}}_2}^{\′}(x\·y))}{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}\underset{y=1}{\overset{m}{\mathrm{\Σ}}}{\left(V_{{\mathrm{\ω}}_2}^{\′}(x\·y)\right)}^2}$ 计算加权因子a，

然后对数据进行加权求差处理，得到边界电压变化ΔV。

$\mathrm{\ΔV}=a\·V_{{\mathrm{\ω}}_2}^{\′}-V_{{\mathrm{\ω}}_1}^{\′}$

$\mathrm{\ΔV}=\left[\begin{array}{c}2.249e-04-3.140e-05i\\ 0.002279-0.0001950i\\ 0.001959-0.0002015i\\ 0.002089-6.580e-05i\\ 0.0009536-2.640e-05i\\ \·\\ \·\\ \·\\ \·\\ \·\\ \·\\ -0.0002328+7.700e-05i\\ -0.0013600+0.0001290i\\ -0.001474+2.830e-05i\\ -0.001338+9.760e-05i\\ 0.001146+1.530e-06i\end{array}\right]$

5、图像重建

由于EIT逆问题求解具有严重的病态性，需要采用合适的正则化方法加以修正，具体选择标准的Tikhonov正则化方法，则Δρ表示为：

$\mathrm{\Δ\ρ}=a\·{\left(\frac{{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}}{{\mathrm{\ρ}}_{\mathrm{\α}}}\right)}^2\·{(S^{T}S+\mathrm{\λR})}^{-1}{S}^T\·\mathrm{\ΔV}$

其中，R为正则化矩阵，λ为正则化系数，本实施例中，R=S^TS,λ=0.1。

$\mathrm{\Δ\ρ}=\left[\begin{array}{c}-0.1084-0.01253i\\ -0.08415-0.01108i\\ -0.05344-0.007547i\\ -0.03495-0.005902i\\ -0.03700-0.007158i\\ \·\\ \·\\ \·\\ \·\\ \·\\ \·\\ 0.0005174-0.0001294i\\ 0.001649+0.0003817i\\ 0.003497+0.0008525i\\ 0.003835+0.0008040i\\ 0.002888+0.0004465i\end{array}\right]$

通过上式即可重建出EIT图像（实部图像和虚部图像），重建结果如图2所示。为了与传统重建算法相比较，本实例中也利用传统的算法重建EIT图像，传统算法的重建结果如图3所示。

从图2和图3中可以看出，加权求差准静态成像算法的重建结果清晰，而传统的图像重建算法的重建结果中存在大量伪影，甚至淹没了目标（虚部图像）。

Claims (6)

1.一种准静态电阻抗成像中的图像重构方法，其特征在于，包括以下步骤：

1）采用电阻抗成像方法，对待测目标分别在频率为和

的电流激励下进行激励，采集两次电流激励下的边界电压

和

然后对所采集的数据进行加权求差处理，得到不同频率激励下的目标区域阻抗变化差异所引起的边界电压变化ΔV：

$\mathrm{\ΔV}=a\·V_{{\mathrm{\ω}}_2}^{\′}-V_{{\mathrm{\ω}}_1}^{\′}$

其中，a为两种不同频率下成像区域内阻抗的变化比例系数，

n为电极数；m为一次激励所产生的边界测量电压数；

2）然后计算不同频率激励下的目标区域阻抗分布变化Δρ：

其中，为激励频率为对应的敏感性矩阵；

利用不同频率激励下的目标区域阻抗分布变化Δρ重建出EIT图像。

2.如权利要求1所述的准静态电阻抗成像中的图像重构方法，其特征在于，根据检测对象估计其场域内阻抗分布为ρ_α，则Δρ为：

$\mathrm{\Δ\ρ}=a\·{\left(\frac{{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}}{{\mathrm{\ρ}}_{\mathrm{\α}}}\right)}^2{S}^{-1}\·\mathrm{\ΔV}$

其中，

和

相等或接近，real(ρ)和imag(ρ)表示ρ的实部和虚部，为实常数，敏感矩阵S的元素通过场域内阻抗分布ρ_α计算：

3.如权利要求2所述的准静态电阻抗成像中的图像重构方法，其特征在于，所述的ρ_α的值与实际检测对象的阻抗分布值成比例，使

为实常数，或是虚部与实部相比较很小，不影响成像结果。

4.如权利要求2所述的准静态电阻抗成像中的图像重构方法，其特征在于，在利用Δρ重建出EIT图像时，还需采用正则化方法加以修正。

5.如权利要求4所述的准静态电阻抗成像中的图像重构方法，其特征在于，采用的Tikhonov正则化方法加以修正：

$\mathrm{\Δ\ρ}=a\·{\left(\frac{{\mathrm{\ρ}}_{{\mathrm{\ω}}_1}}{{\mathrm{\ρ}}_{\mathrm{\α}}}\right)}^2\·{(S^{T}S+\mathrm{\λR})}^{-1}{S}^T\·\mathrm{\ΔV},$ 其中，R为正则化矩阵，λ为正则化系数。

6.如权利要求1所述的准静态电阻抗成像中的图像重构方法，其特征在于，所述的

低于

小于10KHz，大于10KHz。

Images