体内植入电子器件具有光明的、广泛的临床应用前景。在目前的植入电子器件应用中，主要以电池技术和磁感应技术供能。电池供能的最大问题是使用寿命受到限制，电能耗尽后需通过外科手术更换电池；而植入电池占植入器件总体积的50%以上，是妨碍微型化的关键因素^[1-2]。电磁感应技术虽可将体外电能传递到体内且应用广泛，但会对附近的医疗设备产生射频干扰，而其体内感应线圈限制了植入器件微型化^[3]。体导电能量传递技术是植入电子器件能量供给技术的一个重要研究方向。该方式可将丰富的体外电能传递到体内植入器件，是一种有效而稳定的能量供给方案。体导电能量传递系统^[4-6]如图 1所示，主要包括体外电能供给电路、皮肤组织及其两侧的电极组成的皮肤电极单元、体内能量转换电路和功能植入器件，是一个复杂的场路耦合系统。由于皮肤电极单元的电极电流是集总电流，文献[2]提出用X型等效二端口网络(图 2)等效替换皮肤电极单元，建立电路模型，研究体导电能量传递系统。因此，能量传递效率^[7]等特性与该网络的阻抗参数有关。目前，主要通过实验测量的方法获得该网络阻抗参数，尚无从理论计算方法的角度来获取。

笔者研究X型等效二端口网络参数的计算方法，为设计能量传递效率最优的皮肤。

电极单元奠定基础，提出的计算方法如下：

1) 在电极形状和布局一定的条件下，外电极施加正弦电压，通过电磁场原理计算端口电流电压。

2) 通过快速傅里叶变换(FFT)提取端口电流电压的幅值和相位。

3) 在上述激励电压和负载电流确定的条件下，根据X型等效电路(图 2)求取其阻抗参数。

1 场路耦合法计算端口电流电压

为了获得体导电中皮肤电极单元端口处的电压、电流参数，建立皮肤电极单元的场路耦合模型。

1.1 场路耦合物理模型

体导电“场”“路”耦合物理模型如图 3所示，Ω₁、Ω₂、Ω₃、Ω₄为电极区域；Ω₅为皮肤组织区域，含有传导电流和位移电流但不含源电流；S₁、S₂、S₃、S₄为电极与电路的连接边界；S₁₅、S₂₅、S₃₅、S₄₅为电极与皮肤组织的内部交界面；u_in为施加在电极上的交变电压；i₁为注入电极的外电流；R₀为体内电路部分的等效电阻；u_out为其输出电压；i₂为该电路的电流。

体外电压源在电极上施加交变电压，在皮肤组织区域产生交变电场。该电场引起皮肤组织中的带电粒子定向移动，从而在体内、外电极间产生传导电流和位移电流，体内电极接受传递的能量并向体内植入器件供电。

在皮肤电极区域按“场”^[8]的观点进行电磁场计算，而在电路上将这些电极与电源、电容、电阻等集总参数相连。这样以电压源为激励，电流作为未知量进行处理符合体导电能量传递的实际情况。

1.2 准静态电场有限元求解

体导电能量传递系统的交流电压源频率低于100 kHz时，电磁波的波长大于3 000 m，远大于人体的3维尺寸。因此，把皮肤内的时变电磁场简化为准静态电场^[9]。

假设生物组织是各向同性的均匀介质，在皮肤、电极区域内，用磁场强度H和电场强度E表示的场方程和边界条件为

$ \left\{ {\begin{array}{*{20}{l}} {\nabla \times H = \sigma E + \frac{{\partial D}}{{\partial t}},}\\ {\nabla \times E = 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \mathit{\Omega} }\\ {\nabla \cdot D = \rho ,} \end{array}} \right. $

(1)

$ \int_S { - n\left( {J + \frac{{\partial D}}{{\partial t}}} \right)} = {I_1},{S_1}、{S_2}; $

(2)

$ \int_S { - n\left( {J + \frac{{\partial D}}{{\partial t}}} \right)} = {I_2},{S_3},{S_4}; $

(3)

$ n \cdot E = 0,{S_5}; $

(4)

$ n \cdot \left( {{E_1} - {E_2}} \right) = 0\;{S_{15}}、{S_{25}}、{S_{35}}、{S_{45}}。$

(5)

引入标量电位φ来描述电场强度，即

E=-∇φ，则上述场方程可写成：

$ \left. { - \frac{{{\xi _0}{\xi _r}{\nabla ^2}\left( {\partial \varphi } \right.}}{{\partial t}}} \right) - \nabla \cdot \left( {\sigma \nabla \varphi } \right) = 0,\mathit{\Omega }; $

(6)

$ \int_s {n \cdot \left( {\sigma \left( {\nabla \varphi } \right) + {\xi _0}{\xi _r}\nabla \frac{{\partial \varphi }}{{\partial t}}} \right){\rm{d}}s} = {I_1},{S_1},{S_2}; $

(7)

$ \int_s {n \cdot \left( {\sigma \left( {\nabla \varphi } \right) + {\xi _0}{\xi _r}\nabla \frac{{\partial \varphi }}{{\partial t}}} \right){\rm{d}}s} = {I_2},{S_3},{S_4}; $

(8)

$ \frac{{\partial \varphi }}{{\partial n}} = 0,{S_5}; $

(9)

$ \frac{{\partial {\varphi _1}}}{{\partial n}} - \frac{{\partial {\varphi _2}}}{{\partial n}} = 0,\;{S_{15}},{S_{25}},{S_{35}},{S_{45}}。$

(10)

式中：Ω为整个求解区域；σ为电导率；ξ_r为相对介电常数。对上述电场方程(6)~(10)采用加权余量法进行伽辽金有限元离散，并整理，式(6)转换成式(11)：

$ \begin{array}{*{20}{c}} {\int_\mathit{\Omega } {\nabla \cdot \left[ {\left( {\sigma \nabla \varphi } \right){N_j}} \right]{\rm{d}}v} - \int_\mathit{\Omega } {\left[ {\left( {\sigma \nabla \varphi } \right)\nabla \cdot {N_j}} \right]{\rm{d}}v} + }\\ {\int_\mathit{\Omega } {\nabla \cdot \left[ {{N_j}{\xi _0}{\xi _r}\nabla \left( {\frac{{\partial \varphi }}{{\partial t}}} \right)} \right]{\rm{d}}v} - }\\ {\int_\mathit{\Omega } {\left[ {{\xi _0}{\xi _r}\nabla \left( {\frac{{\partial \varphi }}{{\partial t}}} \right)\nabla \cdot {N_j}} \right]{\rm{d}}v} = 0\left( {j = 1,2, \cdots ,NP} \right)。} \end{array} $

(11)

其中：N_j为权函数；NP为节点总数。

根据边界条件式(7)、式(8)、式(9)和式(10)，经过简化，式(11)最后成为

$ \begin{array}{*{20}{c}} {2\int_{S1} {\frac{{{I_1}}}{{{S_1}}}{N_j} \cdot n{\rm{d}}s} + 2\int_{S3} {\frac{{{I_2}}}{{{S_3}}}{N_j} \cdot n{\rm{d}}s} - \int_\mathit{\Omega } {\left( {\sigma \nabla \varphi } \right)\nabla \cdot {N_j}{\rm{d}}v} - }\\ {\int_\mathit{\Omega } {\left( {{\xi _0}{\xi _r}\nabla \frac{{\partial \varphi }}{{\partial t}}} \right)\nabla \cdot {N_j}{\rm{d}}v} = 0。} \end{array} $

(12)

式(12)在整个求解域内空间离散化，可以得到一个时间的一阶微分方程组，将其矩阵形式表示为

$ \mathit{\boldsymbol{K}}\varphi + \mathit{\boldsymbol{M}}\frac{{\partial \varphi }}{{\partial t}} + \mathit{\boldsymbol{CI}} = 0, $

(13)

式中: K为劲度矩阵; M为阻尼矩阵; C、D为耦合劲度矩阵。

1.3 场路耦合电路方程

考虑离子型导电的特点和电极极化效应，体导电能量传递系统采用交流激励源，体外场路耦合电路模型如图 4所示。

体外电路回路方程为

$ {u_{{\rm{in}}}} = i{R_1} + {\varphi _1} - {\varphi _2}, $

(14)

式中:u_in为体外激励电压; i₁为回路电流; R₁为电源等效内阻; $φ$₁、$φ$₂分别为电极1、2的电位。写成矩阵即

$ D\varphi + {R_1}I = {u_{{\rm{in}}}}, $

(15)

体内场路耦合模型如图 5所示。

R₂为等效串联内阻，体内回路电路方程为

$ {i_2}{R_2} = {\varphi _3} - {\varphi _4}, $

(16)

其中:i₂为体内回路电流；$φ$₃、$φ$₄分别为体内电极3、4的点位。

写成矩阵即

$ E\varphi = - {R_2}I = 0。$

(17)

1.4 场路耦合数学模型

将式(13)与式(15)、式(17)联立, 可得体导电的场路耦合模型为

$ \left[ {\begin{array}{*{20}{c}} K&0&C\\ D&{{R_1}}&0\\ E&{ - {R_2}}&{ - 1} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \varphi \\ I\\ 0 \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} M&0&0\\ 0&0&0\\ 0&0&0 \end{array}} \right]\frac{{\rm{d}}}{{{\rm{d}}t}}\left[ {\begin{array}{*{20}{c}} \varphi \\ I\\ 0 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0\\ {{u_{{\rm{in}}}}}\\ 0 \end{array}} \right], $

(18)

将其写成标准形式：

$ \mathit{\boldsymbol{PX}} + \mathit{\boldsymbol{Q}}\frac{{{\rm{d}}\mathit{\boldsymbol{X}}}}{{{\rm{d}}t}} = \mathit{\boldsymbol{f}}, $

(19)

其中：

$ \mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} K&0&C\\ D&{{R_1}}&0\\ E&{ - {R_2}}&{ - 1} \end{array}} \right],\mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} \varphi \\ I\\ 0 \end{array}} \right], $

$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} M&0&0\\ 0&0&0\\ 0&0&0 \end{array}} \right],\mathit{\boldsymbol{f}} = \left[ {\begin{array}{*{20}{c}} 0\\ {{u_{{\rm{in}}}}}\\ 0 \end{array}} \right]; $

式(19)为时域模型，求解时需进行时间离散。笔者采用中心差分方法进行，其离散格式为

$ \frac{{Q\left( {{\varphi _{n + 1}} - {\varphi _n}} \right)}}{{\Delta t}} + \frac{{P\left( {{\varphi _{n + 1}} - {\varphi _n}} \right)}}{2} = \frac{{{f_{n + 1}} + {f_n}}}{2}。$

(20)

将式(19)代入式(20)，整理得：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\frac{M}{{\Delta t}} + \frac{K}{2}}&0&{\frac{C}{2}}\\ {\frac{D}{2}}&{\frac{{{R_1}}}{2}}&0\\ {\frac{E}{2}}&{ - \frac{{{R_2}}}{2}}&{ - \frac{1}{2}} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} \varphi \\ I\\ 0 \end{array}} \right]_{n + 1}} = \\ \left[ {\begin{array}{*{20}{c}} {\frac{M}{{\Delta t}} - \frac{K}{2}}&0&{ - \frac{C}{2}}\\ { - \frac{D}{2}}&{ - \frac{{{R_1}}}{2}}&0\\ {\frac{E}{2}}&{\frac{{{R_2}}}{2}}&{\frac{1}{2}} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} \varphi \\ I\\ 0 \end{array}} \right]_n} + \left[ {\begin{array}{*{20}{c}} 0\\ {\frac{{{u_{in,n + 1}} + {u_{in,n}}}}{2}}\\ 0 \end{array}} \right]。\end{array} $

(21)

式(21)即为体导电能量传递3维瞬态电磁场的场路耦合计算公式。所推导的场路耦合方程虽然增加了几个电流变量，但是方程系数的稀疏对称性的优点却得到了很好的保持；只要知道激励电压条件u_in，便可将此时的标量电位$φ$及电极处的瞬变电流i同时求解出来。

2 FFT提取端口电流电压的幅值和相位

场路耦合模型已仿真得到皮肤电极单元各端口电压、电流的离散数据，为获取它们的幅值和相位，需对其进行谱分析。考虑到全相位FFT(apFFT)谱分析与有限FFT变换相比具有频谱泄漏少的优点^[10-11]，且在无噪声的条件下相位无需要校正^[12]，故选用前者。

现以提取电极电压的幅值和相位为例，apFFT谱分析的具体过程如下：场路耦合模型已得到的电压离散数据为：x(-20), x(-19), …x(0), …, x(20)，共41个数据。离散数据来自体导电的单频信号，故可以列出以x(0)为中心的N个N维向量，有：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_0} = {{\left[ {x\left( 0 \right),x\left( 1 \right), \cdots ,x\left( {N - 1} \right)} \right]}^{\rm{T}}},}\\ {{\mathit{\boldsymbol{x}}_1} = {{\left[ {x\left( { - 1} \right),x\left( 0 \right), \cdots ,x\left( {N - 2} \right)} \right]}^{\rm{T}}},}\\ { \cdots \cdots }\\ {{\mathit{\boldsymbol{x}}_{N - 1}} = {{\left[ {x\left( { - N + 1} \right),x\left( { - N + 2} \right), \cdots ,x\left( 0 \right)} \right]}^{\rm{T}}}。} \end{array} $

其中，N=21。

将每个向量进行循环移位，把样本点x(0)移到首位，则可得到另外的N个N维向量。对准x(0)相加并取其平均值，则可得到全相位数据向量：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_{ap}} = \frac{1}{N}\left[ {Nx\left( 0 \right) + \left( {N - 1} \right)x\left( 1 \right) + x\left( { - N + 1} \right) + \cdots + } \right.}\\ {{{\left. {x\left( {N - 1} \right) + \left( {N - 1} \right)x\left( { - 1} \right)} \right]}^{\rm{T}}}。} \end{array} $

(22)

对式(22)进行FFT变换，即得到apFFT谱：

$ {X_{ap}}\left( k \right) = \frac{A}{{{N^2}}}\frac{{{{\sin }^2}\left[ {{\rm{ \mathsf{ π} }}\left( {\beta - k} \right)} \right]}}{{{{\sin }^2}\left[ {{\rm{ \mathsf{ π} }}\left( {\beta - k} \right)/N} \right]}}{{\rm{e}}^{j{\varphi _0}}},k = 0,1, \cdots ,N - 1。$

其幅值为

$ A = \sum\limits_{k = {k_0} - 2}^{{k_0} + 2} {\left| {{X_{ap}}\left( k \right)} \right|} , $

初相位为

$ {\varphi _0} = {\rm{ang}}\left[ {{X_{ap}}\left( {{k_0}} \right)} \right]。$

按照上述方法同样提取各个电极处电流信号的幅值和初相位，由此可进一步得到电流电压的实部和虚部，为接下来用变载法计算X等效电路阻抗参数提供必要的数据。

3 变载法计算X型等效电路阻抗

目前可以用来计算X等效模型Z参数的方法有开路法、短路法、最优化方法和参数估计法等。其中，开路法和短路法计算体导电网络参数均会改变皮肤电极域的流线分布结构，因此该方法计算的网络参数不能反映体导电系统正常工作时皮肤电极单元域内的结构特征，导致计算结果不准确。另外，该方法还存在多次改变电路结构和多次提取数据的缺点。参数估计法计算体导电网络参数时数学处理较复杂。最优化方法不仅数学处理复杂，且得到的解不唯一。因为皮肤电极单元参数的计算为一非线性多维无约束最优化问题^[13]，约束条件不足，该方法往往能求得多个局部最优解，但难以确定全局最优解。

为克服以上方法的不足，笔者引入了一种可以准确计算X等效模型阻抗参数的方法—基于场路耦合的变载法。变载法通过改变负载端负载值，建立方程组计算阻抗参数。该方法既可以避免改变激励源的接线方式，又可以减少方程的个数，使测量数据的数学处理变得简单。

具体来讲，在1、2端串联接入电压源和电阻R_in，在3、4端接入负载R。为简化分析计算，仅考虑激励源内阻和负载均为纯电阻的情况。由此，图 2变为图 7所示。

当负载R为R₁时，电流电压分别记为I₁⁰、I₂⁰、U₁⁰、U₂⁰；改变负载为R₂时，另记为I₁¹、I₂¹、U₁¹、U₂¹。根据二端口网络知识，得到Z参数方程组(23)和(24)。

$ \left\{ \begin{array}{l} U_1^0 = {Z_{11}}I_1^0 + {Z_1}2I_2^0,\\ U_2^0 = {Z_2}1I_1^0 + {Z_{22}}I_2^0; \end{array} \right. $

(23)

$ \left\{ \begin{array}{l} U_1^1 = {Z_{11}}I_1^1 + {Z_{12}}I_2^1,\\ U_2^1 = {Z_{21}}I_1^1 + {Z_{22}}I_2^1; \end{array} \right. $

(24)

令Z_pq=Z_pqR+jZ_pqI(p, q=1, 2)，将方程组(23)、(24)按实部虚部分解，可得到含8个未知数8个方程的方程组。即

$ \begin{array}{l} \left\{ \begin{array}{l} U_{mR}^n = {Z_{m1R}}I_{1R}^n - {Z_{m1I}}I_{1I}^n - {Z_{m2R}}I_{2R}^n + {Z_{m2I}}I_{2I}^n,\\ U_{mI}^n = {Z_{m1R}}I_{1I}^n - {Z_{m1I}}I_{1R}^n - {Z_{m2R}}I_{2I}^n + {Z_{m2I}}I_{2R}^n; \end{array} \right.\\ \;\;\;\;\left( {其中\;m = 1,2;n = 0,1} \right) \end{array} $

(其中m=1, 2；n=0, 1)

因为皮肤电极单元对称，所以Z₁₂=Z₂₁。于是方程组成为含6个未知数8个方程的超定线性方程组。代入电流电压的实部和虚部值，解此方程组，可以求出皮肤电极单元等效二端口网络阻抗参数Z₁₁、Z₁₂(Z₂₁)和Z₂₂的实部和虚部。再根据以下关系：

$ \left\{ \begin{array}{l} {Z_{11}} = {Z_1} + {Z_2} + {Z_5},\\ {Z_{22}} = {Z_3} + {Z_4} + {Z_5},\\ {Z_1} = {Z_2},\\ {Z_3} = {Z_4},\\ {Z_{12}} = {Z_{21}} = {Z_5}。\end{array} \right. $

即可求出X等效模型的各个阻抗参数。

4 仿真计算和验证 4.1 仿真计算

本仿真实验采用COMSOL Multiphysics3.5a软件搭建实验平台，实现场路耦合的仿真计算。考虑心脏起博器植入胸腔的临床应用，建立皮下植入器件的电极皮肤单元的物理模型，如图 8所示。图中皮肤厚度为6 mm，半径为40 mm，根据意大利应用物理研究院的人体组织的电介特性选择皮肤组织的电导率和相对介电常数，见表 1^[14-15]；电极材料选用导电率为50 000 s/m的AgCl，厚度为1 mm，内外电极紧贴皮肤组织(其距离为皮肤组织厚度)；外激励电压源U_in=20 V，内阻R_in=500 Ω，内电极负载电阻R取R₁=1 000 Ω，R₂=2 000 Ω。因耦合电路是通过spice语言设置，所以图中模型的电路部分没有显示出来。

为了考察电源频率和电极分布对X等效模型各个阻抗参数的影响，进行了下述3类仿真试验。

图 9为电极半径r=5 mm、电极间距d=10 mm、频率取不同值时得到的X等效模型的阻抗参数。可以看出，当频率增大时，阻抗参数的实部和虚部都减小，但当频率达到4 kHz以后，实部变化较缓慢，而虚部变化也平缓了许多。因此，激励频率应选择大于4 kHz。

图 10为频率f=10 kHz、电极半径r=5 mm、电极间距取不同值时得到的X等效模型的阻抗参数。当电极间距增大时，Z₁和Z₃基本保持不变，而Z₅随着间距的增大而增大，与电极间距的增大基本上呈线性变化关系。因此，设计电极时，应适当增大电极间距来增大Z₅，减小横向电流，提高电流传递效率。

图 11为频率f=10 kHz、电极间距d=10 mm、电极半径取不同值时得到的X等效模型的阻抗参数。当半径很小时，阻抗参数很大，随着电极半径的增大，Z₁、Z₃、Z₅都减小；当半径增大到5 mm以后，阻抗参数变化缓慢。因此，宜将电极半径设计在5 mm以上。

4.2 结果验证

采用变载法得到的X型等效电路参数搭建电路，比较该电路和场路耦合在相同激励和负载情况下的端口电流电压，可以检验变载法计算的X型等效电路参数的合理性。图 12为PSIM软件搭建的验证电路，图中电极皮肤阻抗单元的参数是电极间距d=10 mm、电极半径r=5 mm、频率f=10 kHz情况下采用变载法的计算值。激励和输入输出电阻与COMSOL设置保持一致：激励电压源U_in=20 V，频率f=10 kHz，输入电阻R_in=500 Ω，输出电阻R=1 000 Ω。计算的端口电流电压见表 2，相对误差小于0.12%，说明变载法计算的X型等效电路参数是合理的。此外，体内电极电流与体外电极电流之比(即电流传递效率)为58%。

5 结语

笔者通过建立场路耦合模型、运用apFFT谱分析和二端口网络的相关理论，采用变载法成功计算出皮肤电极单元X等效电路模型的阻抗参数，并验证了计算结果的正确性。在此基础上研究了这些参数随供能频率、电极间距和电极半径的变化特性。为进一步从电路角度分析和优化体导电能量传递系统的性能提供了基本参数。