CN104434099A - 利用线性双旋度方程进行磁热声成像的电阻率重建方法 - Google Patents

Abstract

一种利用线性双旋度方程进行磁热声成像的电阻率重建方法，第一步获取断层圆周扫描的导电物体的热声信号；第二步根据磁热声的声压波动方程，利用时间反演算法获取导电物体任意断层面上的热声源分布；第三步结合电流连续性定理，引出矢量电位，对导电物体的电阻率进行初值假定；第四步利用电阻率的初值求解矢量电位空间分量，然后把求解的矢量电位空间分量代入热声源分布，得到更新后的电导率分布；第五步对更新后的电阻率和初值设定的电阻率比对，查看是否满足给定的相对误差，如果满足相对误差则更新后电阻率即为求解的电阻率，否则重新进行第三步的迭代，把更新后的电阻率作为初值电阻率进一步求解更新的电阻率，得到导电物体的电阻率分布。

Description

利用线性双旋度方程进行磁热声成像的电阻率重建方法

技术领域

本发明涉及一种电阻率图像的重建方法，特别涉及一种利用线性双旋度方程进行磁热声成像的电阻率重建方法。

背景技术

目前传统电阻抗成像技术的灵敏度和空间分辨率不高，主要因为电阻抗成像通常采用频率较低的电磁波作为激励，由于波长远远大于成像体，导致电磁场探测对比度高，但分辨率低。毋庸置疑，单一场都有其物理局限性。因此多物理场成像技术受到越来越多的关注，即将一种物理场作用于生物组织，转换为另一种物理场进行检测，由一种物理场提供分辨率，另一种物理场提供对比度，实现对比度和分辨率的同时提高。电磁场和超声相结合的多物理场成像技术正是考虑到电磁场对人体组织电导率的高对比度和超声波探测的高分辨率特性，成为人们的研究热点，磁热声成像作为一种新型的多物理场成像技术在最近一年受到重视。

磁热声成像是由新加坡南洋理工大学在2013年首次提出的新型的电阻抗成像方法，通过对导电物体施加低于20MHz的交变磁场，在导电物体内部产生感应电场，进而产生焦耳热，激发热弹性的声信号，检测声信号进行成像。该方法是一种以交变磁场作为激励源，基于生物组织内部焦耳热吸收率的差异，以超声作为信息载体的无损生物医学影像技术。与微波热声成像技术相比，激励源的频率降低，可以深入到导电体的更深处，使磁热声图像扩展到人体组织的深层。由测量的超声信号到电阻率的重建分为两个过程，首先由测量的超声信号重建热声源分布，然后利用热声源分布重建电阻率分布，目前的相关文献和专利只重建了热声源(S＝ρJ²)，这里E为电场强度的空间分量)，而没有提及电阻率ρ的重建。显然，电场强度E与电阻率ρ的分布有关，从热声源S中重建出电阻率ρ是非常困难的。

发明内容

本发明的目的是克服现有的磁热声成像方法无法给出电阻率分布的不足，提出一种利用线性双旋度方程进行热声源分布重建电阻率分布的重建方法。本发明利用线性双旋度方程进行磁热声成像，可以精确的重建导电物体的电阻率。

本发明基于磁热声成像原理：利用激励线圈对导电物体施加MHz电流激励，在导电物体内产生焦耳热，进而产生超声信号，利用超声换能器接收超声信号，对接收到的热声信号进行热声信号的采集和处理，得到放大滤波后的热声信号后，采用图像重建算法获取导电物体的电阻率图像。

本发明利用线性双旋度方程进行磁热声成像的电阻率重建方法主要包括五个步骤：第一步获取断层圆周扫描的导电物体的热声信号；第二步根据磁热声的声压波动方程，利用时间反演算法获取导电物体任意断层面上的热声源分布；第三步结合电流连续性定理，引出矢量电位，对导电物体的电阻率进行初值假定；第四步结合欧姆定律，利用电阻率的初值求解矢量电位空间分量，然后把求解的矢量电位空间分量代入热声源分布得到更新后的电导率分布；第五步对更新后的电阻率和初值设定的电阻率进行对比，查看是否满足给定的相对误差，如果满足相对误差则更新后电阻率即为求解的电阻率，否则重新进行第三步的迭代，把更新后的电阻率作为初值电阻率进一步进行求解更新的电阻率，求解出导电物体的电阻率分布。实现导电物体的电阻率的重建。

具体步骤描述如下：

第一步：获取断层圆周扫描的导电物体的热声信号：

首先MHz电流激励源通过激励线圈将MHz电流作用到导电物体上，导电物体由于感应电流的作用产生焦耳热，进而产生热声信号，热声信号通过耦合剂耦合到超声换能器内，超声换能器接收到热声信号后通过超声信号处理、采集子系统进行前置放大、滤二级放大处理后，再经过移动窗口取样积分电路对热声信号进一步进行取样积分处理，获取高信噪比的热声信号并传输到上位机进行存储；

第二步，对导电物体的所有断层进行扫描后，获取每个断层的热声信号，图像重建子系统则利用获取的热声信号进行导电物体的电阻率的重建；

已知磁热声成像的声压波动方程：

${\▿}^{2}p(r,t)-\frac{1}{{c_s}^2}\frac{{\∂}^2}{{\∂t}^2}p(r,t)=-\frac{\mathrm{\β}}{C_P}S\left(r\right){\mathrm{\δ}}^{\′}\left(t\right)---\left(1\right)$

其中r为热声源位置坐标，p(r,t)是声压，c_s为热声源在介质中的传播声速，C_P为导电物体的比热容，β为导电物体的热膨胀系数，δ(t)是狄拉克函数，S(r)为热声源分布，t为时间项。

利用时间反演算法，热声源函数的重建公式为：

其中：R＝|r′-r|，

r'是超声换能器的位置坐标，S_d是超声换能器扫描平面，p′是声压对时间的一阶导数，n是r′位置S_d的单位法线矢量。

热声源分布S可以表示为：

S＝ρJ²＝ρJ·J   (3)

其中ρ为导电物体电阻率，J为导电物体内部电流密度分布。

第三步，考虑电流连续性定理引入矢量电位，有：

$J=\▿\×T---\left(4\right)$

其中为哈密顿算符，T为矢量电位空间分量。

已知法拉第电磁感应定律：

$\▿\×E=-B_1---\left(5\right)$

其中E为电场强度，B₁为磁感应强度。

第四步，考虑欧姆定律J＝σE，将式(4)代入式(5)，有：

$\▿\×\mathrm{\ρ}\▿\×T=-B_1---\left(6\right)$

其中，为哈密顿算符，T为矢量电位空间分量。由式(4)代入式(3)，可得：

$\mathrm{\ρ}=\frac{S}{\▿\×T\·\▿\×T}---\left(7\right)$

S为热声源分布。

将式(7)代入式(6)，有：

$\▿\×(\frac{S}{\▿\×T\·\▿\×T}\▿\×T)=-B_1---\left(8\right)$

在成像体的边界处，施加电绝缘边界条件，有：

$n\×{T|}_{\∂\mathrm{\Ω}}=0---\left(9\right)$

其中，n为法向矢量，为导电物体边界定义。

第五步，利用线性双旋度方程的直接迭代方法重建电阻率，步骤如下：

(1)选取导电物体的某一断层面z＝z₀，断层面上的热声源S(x,y,z₀)已经通过时间反演法得出，导电物体上的热声源S(x,y,z)通过断层数据S(x,y,z₀)在z方向的插值得到；

(2)对导电物体进行空间离散，给出电阻率的初值[ρ]⁰；

(3)将[ρ]⁰代入公式(6)，结合公式(9)进行线性有限元求解，重建得到矢量电位空间分量[T]¹；

(4)将热声源S(x,y,z)和矢量电位空间分量[T]¹代入公式(7)，得到更新后的电阻率[ρ]¹；

(5)利用[ρ]¹替换[ρ]⁰，重复根据步骤(3)和(4)，直到电阻率相对误差满足：ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀

即可停止迭代过程。

式中：ε₀为给定的最大相对误差，ε为相对误差。

通过以上过程则可以重建导电物体的电阻率。

附图说明

图1本发明重建方法所涉及的原始热声信号获取装置结构示意图；

图中：1 电流激励源，2 水槽，3 超声换能器扫描控制器，4 激励线圈，5 导电物体，6 超声换能器，7 超声信号处理、采集子系统，8 上位机。

具体实施方式

以下结合附图和具体实施方式进一步说明本发明。

如图1所示，本发明重建方法所涉及的原始热声信号获取装置结构包括激励模块、水槽和检测模块三部分，激励模块包括电流激励源1和激励线圈4，检测模块包括超声换能器6、超声信号处理、采集系统7、超声换能器扫描控制器3和上位机8，其中导电物体5和超声换能器位于水槽2内的耦合剂中，耦合剂可以为去离子水也可以为绝缘油。激励线圈4与导电物体5之间通过空气隔离，导电物体5位于激励线圈4的均匀磁场区域内；所述的电流激励源1电连接激励线圈4，超声换能器6电连接超声信号处理、采集子系统7，数据信号处理、采集子系统7电连接上位机8，超声换能器扫描控制器3和电流激励源1都通过超声信号处理、采集子系统7内部的控制模块进行控制，上位机8实现采集的超声信号的存储和重建图像的计算和显示。

基于所述的热声信号获取装置，磁热声成像的原理描述为：利用激励线圈4对导电物体施加MHz电流激励，在导电物体5内产生焦耳热，进而产生热声信号，利用超声换能器6接收热声信号，对接收到的热声信号进行超声信号的处理和采集，得到放大滤波后的热声信号后，再经过移动窗口取样积分电路对热声信号进一步进行取样积分处理，获取高信噪比的热声信号，采用图像重建算法获取导电物体5的电阻率图像。

所述的一种利用线性双旋度方程进行磁热声成像的电阻率图像重建方法主要包括五个步骤，第一步先获取断层圆周扫描的导电物体5的热声信号；第二步根据磁热声的声压波动方程，利用时间反演算法获取导电物体5任意断层面上的热声源分布；第三步结合电流连续性定理，引出矢量电位，对导电物体5的电阻率进行初值假定；第四步结合欧姆定律，利用电阻率的初值求解矢量电位空间分量，然后把求解的矢量电位空间分量代入热声源分布得到更新后的电导率分布；第五步对更新后的电阻率和初值设定的电阻率进行对比，查看是否满足给定的相对误差，如果满足相对误差则更新后电阻率即为求解的电阻率，否则重新进行第三步的迭代，把更新后的电阻率作为初值电阻率进一步进行求解更新后的电阻率，求解出导电物体的电阻率分布。实现导电物体的电阻率的重建。从电阻率图像重建算法不难看出，重建算法的过程离不开超声信号的获取，超声信号获取的过程描述如下：

如图1所示，所述的一种利用线性双旋度方程进行磁热声成像的电阻率重建方法重建过程描述如下：

第一步：获取断层圆周扫描的导电物体的热声信号：

首先MHz电流激励源1通过激励线圈4将MHz电流作用到导电物体5上，导电物体5由于感应电流的作用产生焦耳热，进而产生热声信号，热声信号通过耦合剂耦合到超声换能器6内，超声换能器6接收到热声信号后通过超声信号处理、采集子系统7进行前置放大、滤波、二级放大处理后，再经过移动窗口取样积分电路对热声信号进一步进行取样积分处理，获取高信噪比的热声信号并传输到上位机8进行存储；

第二步：对导电物体5的所有断层进行扫描后，获取每个断层的热声信号，上位机8则利用获取的热声信号进行导电物体5的电阻率的重建。电阻率重建过程具体描述如下：

已知磁热声成像的声压波动方程：

${\▿}^{2}p(r,t)-\frac{1}{{c_s}^2}\frac{{\∂}^2}{{\∂t}^2}p(r,t)=-\frac{\mathrm{\β}}{C_P}S\left(r\right){\mathrm{\δ}}^{\′}\left(t\right)---\left(1\right)$

其中r为热声源位置坐标，p(r,t)是声压，c_s为热声源在介质中的传播声速，C_P为导电物体5的比热容，β为导电物体5的热膨胀系数，δ(t)是狄拉克函数，S(r)为热声源分布，t为时间项；

根据声压波动方程，利用时间反演法重建的热声源函数表达式为：

其中：R＝|r′-r|，

r'是超声换能器的位置，S_d是超声换能器所在的平面，p′是声压对时间的一阶导数，n是r′位置S_d的单位法线矢量。

而且热声源S表示为：

S＝ρJ²＝ρJ·J   (3)

其中，ρ为导电物体电阻率，J为导电物体内部电流密度分布。

第三步：考虑电流连续性定理引入矢量电位，有：

$J=\▿\×T---\left(4\right)$

其中为哈密顿算符，T为矢量电位空间分量；

已知法拉第电磁感应定律：

$\▿\×E=-B_1---\left(5\right)$

其中，E为电场强度，B₁为磁感应强度；

第四步：考虑欧姆定律J＝σE，将式(4)代入式(5)，有：

$\▿\×\mathrm{\ρ}\▿\×T=-B_1---\left(6\right)$

其中，ρ为电阻率，为哈密顿算符，T为矢量电位空间分量；

由式(4)代入式(3)，得到：

$\mathrm{\ρ}=\frac{S}{\▿\×T\·\▿\×T}---\left(7\right)$

S为热声源分布。

将式(7)代入式(6)，有

$\▿\×(\frac{S}{\▿\×T\·\▿\×T}\▿\×T)=-B_1---\left(8\right)$

在成像体的边界处，施加电绝缘边界条件，有

$n\×{T|}_{\∂\mathrm{\Ω}}=0---\left(9\right)$

其中n为法向矢量，为导电物体边界定义。

第五步：利用线性双旋度方程进行磁热声电阻率重建的步骤如下：

(1)选取导电物体8的某一断层面z＝z₀，断层面上的热声源S(x,y,z₀)已经通过时间反演法得出，导电物体(5)上的热声源S(x,y,z)通过断层数据S(x,y,z₀)在z方向的插值得到；

(2)对导电物体(5)进行空间离散，给出电阻率的初值[ρ]⁰；

(3)将[ρ]⁰代入公式(6)，结合公式(9)进行线性有限元求解，重建得到矢量电位空间分量[T]¹；

(4)将热声源分布S(x,y,z)和矢量电位空间分量[T]¹代入公式(7)，利用电阻率设定的初值和矢量电位分量求解出更新后的电阻率分布[ρ]¹；

(5)判断电阻率相对误差是否满足ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀，其中，ε₀为给定的最大相对误差，满足即可停止迭代过程，[ρ]¹即为最终求解的电阻率。如果不满足ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀，则利用[ρ]¹替换[ρ]⁰，重复根据步骤(3)和(4)，直到电阻率相对误差满足ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀，其中，ε₀为给定的最大相对误差，即可停止迭代过程，[ρ]¹即为最终求解的电阻率。

通过以上过程则可以重建导电物体的电阻率。

最后将重建的导电物体5的电导率图像在上位机8上显示。

Claims (2)

1.一种利用线性双旋度方程进行磁热声成像的电阻率重建方法，其特征在于：所述的利用线性双旋度方程进行磁热声成像的电阻率重建方法包括五个步骤：第一步获取断层圆周扫描的导电物体的热声信号；第二步根据磁热声的声压波动方程，利用时间反演算法获取导电物体任意断层面上的热声源分布；第三步结合电流连续性定理，引出矢量电位，对导电物体的电阻率进行初值假定；第四步结合欧姆定律，利用电阻率的初值求解矢量电位空间分量，然后把求解的矢量电位空间分量代入热声源分布得到更新后的电导率分布；第五步对更新后的电阻率和初值设定的电阻率进行对比，查看是否满足给定的相对误差，如果满足相对误差则更新后电阻率即为求解的电阻率，否则重新进行第三步的迭代，把更新后的电阻率作为初值电阻率进一步进行求解更新的电阻率，求解出导电物体的电阻率分布。实现导电物体的电阻率的重建。

2.按照权利要求1所述的利用线性双旋度方程进行磁热声成像的电阻率重建方法，其特征在于：所述的电阻率重建方法的步骤具体如下：

第一步：获取断层圆周扫描的导电物体的热声信号：

首先MHz电流激励源(1)通过激励线圈(4)将MHz电流作用到导电物体(5)上，导电物体(5)由于感应电流的作用产生焦耳热，进而产生热声信号，热声信号通过耦合剂耦合到超声换能器(6)内，超声换能器(6)接收到热声信号后通过超声信号处理、采集子系统(7)进行前置放大、滤波、二级放大处理后，再经过移动窗口取样积分电路对热声信号进一步进行取样积分处理，获取高信噪比的热声信号并传输到上位机(8)进行存储；

第二步：对导电物体(5)的所有断层进行扫描后，获取每个断层的热声信号，上位机(8)利用获取的热声信号进行导电物体(5)的电阻率的重建；具体的电阻率重建过程描述如下：

已知磁热声成像的声压波动方程：

${\▿}^{2}p(r,t)-\frac{1}{{c_s}^2}\frac{{\∂}^2}{{\∂t}^2}p(r,t)=-\frac{\mathrm{\β}}{C_P}S\left(r\right){\mathrm{\δ}}^{\′}\left(t\right)---\left(1\right)$

其中，r为热声源位置坐标，p(r,t)是声压，c_s为热声源在介质中的传播声速，C_P为导电物体(5)的比热容，β为导电物体(5)的热膨胀系数，δ(t)是狄拉克函数，S(r)为热声源分布，t为时间项；

根据声压波动方程，利用时间反演法重建的热声源函数表达式为：

其中：R＝|r′-r|，

r'是超声换能器的位置，S_d是超声换能器扫描平面，p′是声压对时间的一阶导数，n是r′位置S_d的单位法线矢量；

而且热声源S表示为：

S＝ρJ²＝ρJ·J    (3)

其中，ρ为导电物体(5)电阻率，J为导电物体(5)内部电流密度分布；

第三步：考虑电流连续性定理引入矢量电位，有：

$J=\▿\×T---\left(4\right)$

其中，为哈密顿算符，T为矢量电位空间分量；

已知法拉第电磁感应定律：

$\▿\×E=-B_1---\left(5\right)$

其中，E为电场强度，B₁为磁感应强度；

第四步：考虑欧姆定律J＝σE，将式(4)代入式(5)，有：

$\▿\×\mathrm{\ρ}\▿\×T=-B_1---\left(6\right)$

其中，ρ为电阻率，为哈密顿算符，T为矢量电位空间分量；

由式(4)代入式(3)，得到：

$\mathrm{\ρ}=\frac{S}{\▿\×T\·\▿\×T}---\left(7\right)$

S为热声源分布；

将式(7)代入式(6)，有：

$\▿\×(\frac{S}{\▿\×T\·\▿\×T}\▿\×T)=-B_1---\left(8\right)$

在成像体的边界处，施加电绝缘边界条件，有

$n\×T{|}_{\∂\mathrm{\Ω}}=0---\left(9\right)$

其中，n为法向矢量，为导电物体边界定义；

第五步：利用线性双旋度方程进行磁热声电阻率重建，步骤如下：

(1)选取导电物体(5)的某一断层面z＝z₀，断层面上的热声源S(x,y,z₀)通过时间反演法得出，导电物体(5)上的热声源S(x,y,z)通过断层数据S(x,y,z₀)在z方向的插值得到；

(2)对导电物体(5)进行空间离散，给出电阻率的初值[ρ]⁰；

(3)将[ρ]⁰代入公式(6)，结合公式(9)进行线性有限元求解，重建得到矢量电位空间分量[T]¹；

(4)将热声源分布S(x,y,z)和矢量电位空间分量[T]¹代入公式(7)，利用电阻率设定的初值和矢量电位分量求解出更新后的电阻率分布^[ρ]¹；

(5)判断电阻率相对误差是否满足ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀，其中，ε₀为给定的最大相对误差，满足即可停止迭代过程，[ρ]¹即为最终求解的电阻率，如果不满足ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀，则利用[ρ]¹替换[ρ]⁰，重复根据步骤(3)和步骤(4)，直到电阻率相对误差满足ε＝||([ρ]¹-[ρ]⁰)/[ρ]⁰||₂≤ε₀，其中，ε₀为给定的最大相对误差，即可停止迭代过程，[ρ]¹即为最终求解的电阻率；

通过以上过程则可以重建导电物体(5)的电阻率；

最后将重建的导电物体(5)的电导率图像在上位机(8)上显示。