组织电阻抗(电导率)对反映病理和生理变化具有重要的临床价值，目前有多个国家竞相开展了电阻抗成像研究。国内外学者主要尝试通过两种途径发展电阻抗成像：1)利用其功能分辨率的临床潜在优势将电阻抗成像推向实用，研究利用电阻抗进行监护的方法^[1]；2)研究与超声、核磁共振等相结合的混合成像方法，解决空间分辨率不高的问题^[2-3]。磁声成像方法(MAT)^[4-6]和磁声电成像(MAET)^[7-9]是两种典型的电阻抗成像和超声相结合的成像方法，融合了电阻抗成像技术和超声扫描成像技术，兼具无创、对比度好、灵敏度高以及空间分辨率高等优点。MAT和MAET两种成像方法互为反模式，MAT利用注入电流(或感应电流)和静磁场激励生物组织产生洛伦兹力振动形成超声波，提取超声信号重建成像体的电导率分布；MAET利用超声波和静磁场产生洛伦兹力，激励生物组织检测电信号重构电阻抗图像。从电工学角度可见，MAT依据的是电动机原理，而MAET依据的则是发电机原理。相比MAET，MAT在国际上是前沿和热点研究方向，研究人员在数学物理模型^[10-11]、声电激励与信号转换实验^[12]、信息检测与图象重建方法^[13]等方面开展了相应的研究工作。

MAET可以看成MAT的反模式，该方法起步较晚，它的具体原理是将一束超声波注入成像体，成像体中的局部离子随超声波的传播而振动，振动的离子在静磁场作用下受到洛伦兹力作用而引起电荷分离，在成像体内形成局部电场，通过贴放在成像体上的接收电极或与成像体非接触的接收线圈检测电信号重构电阻抗图像。MAET作为一种新型的电阻抗成像与超声成像相结合的技术，与传统的电阻抗成像技术相比，具有4方面的优点。

1) 激励方面。磁声电成像采用超声探头遍历置于静磁场中的被测目标体。不同于传统电阻抗的整区激励方式，借助互易定理，利用超声离子振动速度和静磁场信息，分区获取电流密度；而传统的成像方法则是在重建过程中才对目标体进行区域分解，利用的是整体电流密度信息。

2) 测量方面。理论上将仅采用一对接收电极即可实现空间分布区域的电导率重建，这是传统电阻抗成像方法无法做到的。

3) 图像重建方面。磁声电成像方法理论上可以采用超声声速编码技术，精确获取目标体的空间信息，采用超声聚焦技术可以有效地控制图像横向分辨率。

4) 多场耦合方面。磁声电成像采用超声激励与磁场耦合产生的电场作为激励源，结合了电磁探测对比度高和超声成像分辨率高的优点，便于与传统超声成像技术联合使用，有利于信息互补，做出更准确的诊断。

综上所述，磁声电成像是一种有望获得高质量图像的生物电磁成像方法。目前国内外对磁声电的研究还处于起步阶段，尚没有完善的理论模型和系统实验。本文在基于互易定理的理论研究基础上，实现二维模型电流密度和电导率分布的重建图像，同时对超声探头焦斑大小和扫描步长对重建图像质量的影响进行分析。

1 电磁场正逆问题分析 1.1 电磁场正问题研究

设电极在成像体a、b处测得的MAET电压差为${U_{ab}}$，在成像体a、b处外加单位电流后，其内部的电流密度分布为${\mathit{\boldsymbol{J}}_{ab}}$，离子振动速度为${\mathit{\boldsymbol{v}}}$，则根据互易定理，${U_{ab}}$可以表示为^[10]：

$ {U_{ab}} = \int_V {{\mathit{\boldsymbol{J}}_{ab}} \cdot (\mathit{\boldsymbol{v}} \times {\mathit{\boldsymbol{B}}_0}){\rm{d}}V} $

(1)

利用式(1)获得了测量电极a、b处检测到的MAET表面电压信号的定量求解公式。

如式(1)所示，MAET电磁场正问题可以分为两组：1)已知局部电场源为$\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{'}}_{\rm{1}}}$($\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{'}}_{\rm{1}}} = \mathit{\boldsymbol{v}} \times {\mathit{\boldsymbol{B}}_0}$)，求电极a、b两端的表面电压${U_{ab}}$。2)已知电极a、b端注入的单位电流，求解成像体中的电流密度分布${\mathit{\boldsymbol{J}}_{ab}}$。

第一组正问题的几何模型如图 1所示，为厚度可忽略的二维成像体，置于静磁场${\mathit{\boldsymbol{B}}_0}$中，超声脉冲垂直成像体平面入射，忽略成像体的厚度，则认为质点振动速度$\mathit{\boldsymbol{v}}$不变。超声脉冲和静磁场共同作用区为${\mathit{\Omega} _1}$。

局部电场源激发的电场为库仑电场，感应电场非常小可忽略，采用电准静态近似，并引入标量电位$u$，由电流连续性定理，可以导出局部电场源内部${\mathit{\Omega} _1}$区域电场所满足的方程为：

$ \nabla \cdot (\sigma \nabla u) = \nabla \cdot (\sigma \mathit{\boldsymbol{v}} \times {\mathit{\boldsymbol{B}}_0}) $

(2)

在局部电场源外的区域中(记为${\mathit{\Omega} _2}$)有：

$ \nabla \cdot (\sigma \nabla u) = 0 $

(3)

边界条件为：

$ {\left. {\sigma \frac{{\partial u}}{{\partial n}}} \right|_{\partial \mathit{\Omega} }} = 0 $

(4)

由式(2)~式(4)通过数值求解，可以计算出成像体内的电压分布，进一步可计算出电极a、b的电势差${U_{ab}}$。第二组电磁场正问题是把a、b处的测量电极作激励电极，注入单位电流，同理可以导出整个成像体内所满足的方程为：

$ \nabla \cdot (\sigma \nabla u) = \delta (\mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{r}}_a}) - \delta (\mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{r}}_b}) $

(5)

式中，$\mathit{\boldsymbol{r}}$、${\mathit{\boldsymbol{r}}_a}$和${\mathit{\boldsymbol{r}}_b}$分别为场点、电极$a$以及电极$b$的坐标。

边界条件为：

$ \sigma {\left. {\frac{{\partial u}}{{\partial n}}} \right|_{\partial \mathit{\Omega} }} = 0 $

(6)

式(5)和式(6)即为单位电流注入成像体所满足的电磁场边值问题。

1.2 电磁场逆问题研究

由式(1)可以看出，MAET重建可以分为两个步骤：1)重建电流密度${J_{ab}}$，即单位电流通入成像体中内部电流密度分布的重建；2)根据重建的电流密度${J_{ab}}$重建电导率。本文考虑2D样本，采用图 1所示的超声脉冲激励方式，假设超声脉冲激励方向为$z$方向，被测样本位于$x - y$平面内，超声脉冲垂直入射被测样本，不考虑被测样本的厚度，当静磁场${\mathit{\boldsymbol{B}}_0}$的方向为$y$方向时，由式(1)可以把体积分直接转化成面积分，获得每一次超声脉冲和静磁场共同作用区${\mathit{\Omega} _{{\rm{ }}1}}$内，在$x$方向上的平均电流密度为：

$ {\left. {{{\bar J}_{ab}}} \right|_x} = \frac{{{U_{ab}}}}{{\bar v{B_0}ds}} $

(7)

式中，$\bar v$为超声脉冲激励在超声波传播方向上超声脉冲和静磁场共同作用区${\mathit{\Omega} _{{\rm{ }}1}}$内的平均离子振动速度；${B_0}$为静磁场强度的模值；$ds$为不考虑样本厚度情况下，超声脉冲和静磁场共同作用区${\mathit{\Omega} _{{\rm{ }}1}}$的面积。然后其他条件不变，只旋转静磁场，使静磁场的方向为$x$方向，即可获得每一次超声脉冲和静磁场共同作用区${\mathit{\Omega} _{{\rm{ }}1}}$内在$ - y$方向上的平均电流密度${\left. {{{\bar J}_{ab}}} \right|_y}$。结合扫描系统，使超声探头扫描遍历整个成像体，在分别获得所有的超声脉冲和静磁场共同作用区${\mathit{\Omega} _1}$的$x$方向和$ - y$方向的平均面电流密度后，参照磁共振电阻抗成像中的电导率重建算法，采用迭代算法就可以重建电导率分布^[14]。迭代算法的具体过程首先是建立计算电流密度和测量电流密度之间的目标函数，利用有限元进行正问题求解，根据最小平方误差原则极小化目标函数，最后通过迭代获得成像区域内部电导率的分布。目标函数可写为：

$ \begin{array}{l} Q(\mathit{\boldsymbol{\sigma }}) = \sum\limits_{i = 1}^n {\int_{{\mathit{\Omega} _\iota }} {{{\left| {{\mathit{\boldsymbol{J}}^M}(\mathit{\boldsymbol{r}}) - \mathit{\boldsymbol{J}}(\mathit{\boldsymbol{r}})} \right|}^2}} } {\rm{d}}\mathit{\Omega} \\ = \sum\limits_{i = 1}^n {\int_{{\mathit{\Omega} _\iota }} {{{\left| {{\mathit{\boldsymbol{J}}^M}(\mathit{\boldsymbol{r}}) + {\sigma _i}\nabla \mathit{\boldsymbol{U}}} \right|}^2}} {\rm{d}}\mathit{\Omega} } \end{array} $

(8)

式中，$n$为有限元剖分单元数；$\mathit{\boldsymbol{U}}$和$\mathit{\boldsymbol{J}}$分别为通过边值问题获得的电位和电流密度；${\mathit{\boldsymbol{J}}^M}$为由互易定理结合式(7)获得的单位电流注入被测样本后，在超声波探头作用区域获得电流密度。目标函数进一步可以写成：

$ \begin{array}{l} Q(\mathit{\boldsymbol{\sigma }}) = \sum\limits_{i = 1}^n {\int_{{\mathit{\Omega} _i}} {({\mathit{\boldsymbol{J}}^M}(\mathit{\boldsymbol{r}}) + {\sigma _i}\nabla \mathit{\boldsymbol{U}}) \cdot ({\mathit{\boldsymbol{J}}^M}(\mathit{\boldsymbol{r}}) + {\sigma _i}\nabla \mathit{\boldsymbol{U}})} } {\rm{ d}}\mathit{\mathit{\Omega} } \\ = \sum\limits_{i = 1}^n {\int_{{\mathit{\Omega} _i}} {({\mathit{\boldsymbol{J}}^M} \cdot {\mathit{\boldsymbol{J}}^M} + 2{\sigma _i}{\mathit{\boldsymbol{J}}^M} \cdot \nabla \mathit{\boldsymbol{U}} + \sigma _i^2\nabla \mathit{\boldsymbol{U}} \cdot \nabla \mathit{\boldsymbol{U}})} } {\rm{d}}\mathit{\mathit{\Omega} } \end{array} $

(9)

采用最小二乘法对式(9)的目标函数极小化，有：

$ \begin{array}{l} \frac{{\partial Q(\sigma )}}{{\partial {\sigma _i}}} = 2\int_{{\mathit{\Omega} _i}} {({\mathit{\boldsymbol{J}}^M} \cdot \nabla \mathit{\boldsymbol{U}} + {\sigma _i}\nabla \mathit{\boldsymbol{U}} \cdot \nabla \mathit{\boldsymbol{U}})} {\rm{d}}\mathit{\Omega} = 0\\ i = 1,2, \cdots ,n \end{array} $

(10)

因此有：

$ {\sigma _i} = \frac{{ - \int_{{\mathit{\Omega} _i}} {{\mathit{\boldsymbol{J}}^M} \cdot \nabla \mathit{\boldsymbol{U}}} {\rm{ d}}\mathit{\Omega} }}{{\int_{{\mathit{\Omega} _i}} {\nabla \mathit{\boldsymbol{U}} \cdot \nabla \mathit{\boldsymbol{U}}} {\rm{ }}2{\rm{d}}\mathit{\Omega} }}\;\;\;\;i = 1,2, \cdots ,n $

(11)

给定样品初始电导率，通过重复迭代，便可以达到预先给定的精度值，求出电导率分布。

2 电流密度和电导率重建

首先验证电流密度重建和电导率重建的可行性，设计铜片扫描实验，铜片几何尺寸如图 2所示，尺寸为6 cm×10 cm×1 mm，中间两侧切除的部分尺寸为2 cm×3 cm。

通过图 2所述模型，根据式(7)可以计算出各扫描区域内电流密度分布的$x$分量和$y$分量的平均值，将反演计算所得的电流密度平均值矩阵赋值给各几何区域内的几何中心，使用二次插值，插值后的成像矩阵为110×160，然后根据式(11)求出电导率的分布，将各点的电导率值转化为灰度显示，如图 3所示。可以看出，重建的电导率图像能够近似反映实验样本情况。实验平台的搭建及实验结果参见参考文献[9]。

在验证了重建算法的基础上，建立如图 4所示的几何模型，其几何尺寸如图 4所示。设背景区域$\mathit{\Omega} $、圆形几何区域${\mathit{\Omega} _{{\rm{ }}1}}$、方形区域${\mathit{\Omega} _{{\rm{ }}2}}$和椭圆几何区域${\mathit{\Omega} _{{\rm{ }}3}}$的电导率分别为$0.5$、$1$、$0.01$、$0.4{\rm{ }}\rm{S/m}$。

首先设扫描步径为$0.5{\rm{ }}\rm{mm}$，扫描超声脉冲的焦斑大小为$1.6{\rm{ }}\rm{mm}$直径的圆形，不考虑样本厚度，当$\mathit{\boldsymbol{v}} \times {\mathit{\boldsymbol{B}}_0}$为$x$方向时，由${U_{ab}}$通过式(7)可求得在每一个扫描焦斑区域内，电流密度分布$x$分量平均值。旋转静磁场，$\mathit{\boldsymbol{v}} \times {\mathit{\boldsymbol{B}}_0}$为$y$方向时，同理可求出每一个扫描区域内电流密度分布$y$分量平均值平均值，进而可以求出$x\rm{ - }y$面上的电流密度分布，然后对求出的$x\rm{ - }y$面上的所有$x$分量和$y$分量的电流密度进行双线性插值，最后根据式(11)可以求出整个几何模型的电导率的分布，如图 5所示。

图 5a、图 5b和图 5c所示分别为重建的$x$分量电流密度分布图、重建的$y$分量电流密度分布图和重建的电导率分布图。从图中可以看出，对于成像体的厚度忽略的2D样本，采用本文的电流密度重建算法和电导率重建算法，可以很清晰地获取2D被测样本的电流密度分布图和电导率分布图；重建的$x$分量的电流密度分布能够清晰反映出$x$方向的边界，而重建的$y$分量电流密度分布对$y$方向的边界更为敏感。由此得出结论，电流密度的不同方向分量对不同方向的图像横向分辨力有贡献。

3 重建图像质量影响因素分析

为了论证焦斑大小对图像分辨力的影响，仍采用图 4所示模型，扫描步径设为$0.5{\rm{ }}\rm{mm}$，焦斑直径分别为$1.6$、$2$、$4$、$6{\rm{ }}\rm{mm}$，重建出来的电流密度分布在圆形区域的细节放大图如图 6所示。由图中可看出，随着焦斑变大，图像变模糊，边界清晰度下降，即分辨力降低。理论上，磁声电成像的横向分辨率为焦斑大小，因此，该仿真结果与理论吻合，同时直观地说明了随着焦斑变大，横向分辨力降低。脉冲空间宽度是指超声脉冲在介质中传播时，超声波在纵向上的跨度，当脉冲宽度为$t$，介质中声波的传播速度的大小为$v$，则脉冲空间宽度为$vt$，对应图像的轴向分辨率，在2D成像中没有体现。

不改变焦斑大小，改变扫描步径，考查扫描步径对重建图像质量的影响。设焦斑直径为$4{\rm{ }}\rm{mm}$，扫描步径分别为$2$、$1$、$0.5$、$0.1{\rm{ }}\rm{mm}$，对应的重建结果细节图如图 7所示。由图中可看出，扫描步径越小，即扫描点越多，成像越清晰，即扫描步径长度影响图像质量。当焦斑直径为$4{\rm{ }}\rm{mm}$时，扫描步径由$2{\rm{ }}\rm{mm}$减为$1{\rm{ }}\rm{mm}$，重建效应改善明显；但是当扫描步径由$0.5{\rm{ }}\rm{mm}$减为$0.1{\rm{ }}\rm{mm}$，图像没有明显的改善，说明此时已达到该实验的极限分辨力。

4 结论

在基于互易定理的理论研究基础上，本文实现二维模型电流密度分布的重建图像和电导率的图像重建。同时对超声探头焦斑大小和扫描步长对重建图像质量的影响进行定性的分析。综合以上的仿真结果，重建的电流密度图像和电导率图像较好地显示了被测样本形状，区别出不同的电导率。由以上的仿真实验可以得出：1)通过本文的仿真方法，重建图像理论上的横向分辨率由超声探头的焦斑大小决定，焦斑越小横向分辨率越高；纵向分辨率也可以成为轴向分辨率由超声脉冲宽度和超声波在样本中的声波传输速度的乘积决定。2)扫描步径减小，即扫描次数增多，有助于改善成像质量；但当扫描步径小到一定程度后，对图像质量的影响不再明显。