欧洲心脏病学指南提出，中心动脉压(central aortic pulse wave,CAP)作为直接反映心脏和脑血管压力的参数，可能成为终点心血管事件的预测指标^[1].近年来一些前瞻性的临床随访研究证实，CAP与心血管系统状态和特性(如心脏负荷、心肌收缩力、动脉顺应性、外周阻力等)密切相关，这对于预测心血管事件的发生意义重大^[2].

目前CAP可通过有创或无创手段测量.有创测压法不仅要使病人承受很大的创伤和风险，而且要求实施医院有相当好的技术和设备，临床应用存在局限性.无创测量法主要采用传递函数法(transfer function,TF)，它利用外周动脉(如桡动脉、肱动脉)测得的脉搏波经函数转化推测CAP波形与压力值，该方法的准确性和可靠性已由大量实验证实，但其是否受早搏、心房颤动、药物和心率 等因素的影响还有待深入研究^[3].传递函数法通过测量外周动脉的脉搏波(peripheral artery pulse wave,PAP)，利用函数转换关系间接推导出CAP波形与压力值.目前已有类似产品应用于临床^[4].

传递函数(TF)方法在经历了由GTF(generalized transfer function)到 ITF(individualized transfer function)的多次改进后，其精度已有显著提高，但要应用于临床还有一些理论上与实际操作上的问题需要解决.首先，TF与脉搏波传导特性并无直接联系，其中各参数生理意义尚不明确；其次，这一方法仅通过单个部位PAP和一个先验的TF来估计CAP，忽略了心血管系统的经时渐变特性与个体间的差异性；最后，TF的求解不仅需要大量实验人群，还要采用有创测量手段，时间与物质成本较高.近年来，随着盲辨识技术在生物医学领域的广泛应用，J.O. Hahn等^[5]提出了由多路PAP估计CAP的多通道盲辨识法(multi-channel blind system identification,MBSI)，本文将对其进行深入分析与改进.

1 多通道盲辨识算法

MBSI是盲信号处理技术的一个分支，MBSI算法在CAP估计中的应用大致思路是结合体表两点或多点处的PAP信息来估计CAP^[6].该方法将心血管系统视为一个单输入多输出(single-input multi-output,SIMO)系统，如图 1所示.心室射血的压力或流量波形可作为系统的输入，体表多个部位测得的脉搏压力波形可视为系统的多路输出，脉搏波在外周动脉中传导的过程由系统传递函数描述；这样，CAP重建被转换为由多路输出信号估计各通道系数与输入信号问题，即系统辨识领域的MBSI问题^[6].

图 1(Figure 1)

Fig. 1

图1 心血管系统的SIMO模型 Fig. 1 SIMO model of the cardiovascular system

可辨识系统需满足如下三个必要条件^{[5, 7]}：①所有通道间没有公共零极点，即各通道特性要有足够的不同；②输入信号要足够复杂以激励各个通道，其模式个数不小于L+2(L为系统传递函数阶次)，例如，输入序列不能为零、常数或单个正弦信号；③输出信号足够长，每一通道输出的信号长度N>2L+1.

综上所述，MBSI算法克服了传递函数法需建立先验函数或模型的缺点，可实现CAP的无创、个体化测量.Hahn提出的“灰箱”模型还能反映脉搏波的传导特性，可间接应用于心血管系统重要参数的推导^[5].

本文在前人研究的基础上，尝试用MBSI算法估计CAP波形.CAP最直接地反映了心血管系统状态，还可间接用于估计心输出量(cardiac output,CO)、外周阻力(total peripheral resistance,TPR)、主动脉与外周动脉顺应性等心血管系统功能参数^{[5, 8]}.由于FIR函数结构简单且能反映心血管系统特性，而两路PAP估计CAP运算复杂度低，临床可操作性强，因此本文基于心血管系统“T-tube”模型建立双通道“黑箱”模型，研究心血管系统阶次估计与可辨识性等问题.

心血管系统可被视为一多通道动态系统.Hahn提出用非对称“T-tube”模型^[9](见图 2)可仿真人体不同生理病理状态下的脉搏波，该模型结构简单，原理明确，现已应用于心血管系统参数计算、脉搏波反射特性分析等^[5].本文采用这一模型求解脉搏波传导的传递函数，进一步分析CAP与PAP间的关系.

图 2(Figure 2)

Fig. 2

图2 “T-tube”模型结构 Fig. 2 Structure of the‘T-tube’model

“T-tube”模型结构是传输线理论与弹性腔模型结合的产物，主要由传输线与末端阻抗两部分组成，传输线部分表征近心端大血管特性，末端三元弹性腔模型用于模拟外周小动脉的阻力特性^[10].图 2中Γ_P为反射系数，Z_c表示血管的特性阻抗，C为血管终端动脉顺应性，R表示血管外周阻力，R₀的引入可以更好地描述系统高频特性.设主动脉处脉搏波P_A传导至外周动脉形成P_P，其传导时间为 n_A-P，则在S域中P_A与P_P的关系可由式(1)描述.

由z=e^s/F_s,N_A-P=n_A-P · F_s可将式(1)变换到离散域，得到式(2).

Campbell等通过动物实验证明了非对称“T-tube”模型两个通路参数间模型参数随心血管系统生理状态的变化而变化，当血管处于舒张状态时，动脉顺应性增加，外周阻力减小，脉搏波传导速度降低使延迟时间增加；而血管处于收缩状态时，情况正好相反.Campbell等的实验中还测得了狗在多种状态下的脉搏波，将主动脉脉搏波作为模型输入，验证了非对称“T-tube”模型的有效性；在验证过程中，需不断对“T-tube”模型参数进行调整，使输出波形与两路外周动脉实测(肱动脉/颈动脉、股动脉)脉搏波达到最佳拟合^[9].图 3显示了按照“T-tube”模型在三种生理状态下重建出的脉搏波仿真结果，由图可见，血管收缩时血压升高，血管舒张时血压下降，与真实情况相符，也为本文中MBSI算法验证提供了依据.

图 3(Figure 3)

Fig. 3

图3 不同生理状态下“T-tube”模型仿真的脉搏波 Fig. 3 Simulated pulse wave under different physiological states using the‘T-tube’model (a)—正常生理状态； (b)—血管收缩状态； (c)—血管舒张状态.

非对称“T-tube”模型用简单电路元件对心血管系统结构做出了假设，很明显为一灰箱模型.式(2)给出了该模型传递函数G(z,θ)，为一结构基本确定的IIR函数，其未知参数N_A-P由脉搏波传导时间与采样率决定，而a, b则与末端弹性腔模型的参数有关.在 “T-tube”模型基础上，Hahn最先设计出了心血管系统“灰箱”模型^[5]，其结构如图 4所示.

图 4(Figure 4)

Fig. 4

图4 心血管系统“灰箱”模型结构 Fig. 4 Structure of the‘gray-box’model of cardiovascular system

G₁(s),G₂(s)分别为通路1(头通路)与通路2(体通路)的传递函数，式(3)给出了G₁(s),G₂(s)的计算公式，其离散域表达式与式(2)相同，其未知系数可通过MBSI算法由头通路脉搏波P₁(t)与体通路脉搏波P₂(t)确定.

式(3)中G₁^-1(s)，G₂^-1(s)分别为G₁(s)，G₂(s)的逆函数，用于估计主动脉脉搏波P_A(t).MBSI的“灰箱”方法不仅能估计主动脉脉搏波，还可用于计算外周阻力、动脉顺应性、脉搏波传导时间等参数^{[8, 10]}.但该算法难点是IIR分子、分母系数与阶次的确定，利用两路输出仅能辨识出IIR函数分子与分母的乘积D₂(z^-1) · N₁(z^-1)和D₁(z^-1) · N₂(z^-1)，无法确定系统零极点；增加额外通道虽有助于零极点的辨识，但会增加运算与操作复杂度，降低临床实用性^[6].

与“灰箱”模型相比，“黑箱”模型无需心血管系统建模理论支持，仅将各通路传递函数视为结构与系数未知的IIR或FIR函数.Swamy等^[8]认为脉搏波传导特性可满足FIR滤波器的三个假设条件：首先，心血管系统状态在短时间(60~120s)内是稳定的，满足线性时不变特性；其次，由于不同外周动脉处脉搏波形态差异较大，可认为各通道传递函数是互质的，即没有公共零点；最后，假设主动脉脉搏波的高频成分高于滤波器阶次，则FIR滤波器不会将高频噪声引入估计的源信号.Swamy等用FIR函数来描述“黑箱”模型各通道特性，并用动物实验证明了该方法的有效性，主动脉脉搏波估计的整体误差仅为4.8 % .该方法最大的优点是模型结构简单，运算量小^[8].但目前尚无精确算法与准则可用于FIR阶次的估计，辨识过程中需要不断调整阶次.

IIR与FIR函数各有优缺点，在实际应用中需根据系统特性从多方面综合考虑选择.从性能上来说，IIR滤波器的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性；但这个高效率是以相位的非线性为代价的，选择性越好，则相位非线性越严重.相反，FIR滤波器可以得到严格的线性相位，然而由于FIR滤波器的极点固定在原点，所以只能用较高的阶数达到高的选择性.对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高5~10倍，运算过程复杂，信号延时也大.当对选择性和线性的要求同等重要时，IIR滤波器就必须加全通网络进行相位较正，同样会增加滤波器的复杂性.

从结构上比较，IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定；另外，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引入寄生振荡.相反，FIR滤波器主要采用非递归结构，不论在理论上还是实际运算中都不存在稳定性问题，运算误差也较小.此外，采用快速傅立叶变换法求解FIR滤波器，在阶数相同时，可大大提高运算速度.

整体看来，IIR与FIR滤波器作为系统传递函数时可达到同样的滤波效果；FIR主要的优势是其线性相位、简单性与稳定性，而IIR的长处在于可用低阶非线性相位设计实现频域滤波.在心血管系统特性描述上，基于“T-tube”模型的IIR传递函数有着明显优势，为满足相同的幅频响应，FIR和IIR所需的阶数相差悬殊.但就MBSI算法而言，FIR系数的盲辨识更容易实现.权衡这两大条件，本文尝试用FIR函数来描述脉搏波传导特性，仿真数据选用图 3中的三例典型脉搏波，具体步骤如下.

1) 频域求解两通路传递函数.

两通路传递函数可通过系统输入与输出信号间的互功率谱密度和输入信号的功率谱密度的比值由Welch’s平均周期法估计得到.其中，窗函数采用汉明窗.将主动脉脉搏波P_A作为输入，两路外周脉搏波P₁与P₂分别作为输出，可求得FIR类型传递函数TF_A1与TF_A2.

2) 两通路FIR系数求解.

两通路FIR系数可由频域传递函数的幅值与相位信息求解时域FIR模型得到.设n与m分别为IIR或FIR函数分子与分母阶次；输出b和a即为IIR或FIR的分子与分母系数.将步骤1中求得的TF_A1,TF_A2分别代入，设置m=1,n=5,6,…,40，采用最小二乘方法可得阶次分别为5到40 的两路FIR传递函数各36组，系数b_ij中i=1,2表示通道的选择，j=5,6,…,40为FIR阶次.

3) 还原两路输出波形，评价FIR模型可行性.

将FIR的系数b_ij代入FIR差分数字滤波器后得到此FIR滤波器的输出P_ij.各阶次FIR性能可由系统输出P_1j,P_2j与原始波形P₁,P₂间的误差来评价.评价标准波形畸变率(waveform distortion,WFD)与整体能量差(total energy difference,TED)的定义如下：

图 5显示，两波形间整体误差很小，尤其是收缩压与舒张压误差均在3mmHg以内，已达临床可接受范围.“T-tube”模型对应为IIR数字系统.

图 5(Figure 5)

Fig. 5

图5 “T-tube”输出波形(IIR输出)与FIR估计波形比较 Fig. 5 Comparison of the output waves from ‘T-tube’ model and FIR (a)—正常生理状态； (b)—血管收缩状态； (c)—血管舒张状态.

图 6显示了步骤3的误差分析结果，其中Pulse1,Pulse2,Pulse3分别对应图 3中正常生理状态、血管收缩状态以及血管舒张状态的脉搏波.随FIR阶次n的增加，WFD与TED均逐渐减小；但WFD不受脉搏波均值影响，更能反映波形间形态上的差异.由图中虚线可知，当n在15~20范围内时误差下降幅度开始减小；而n>20时，三种生理状态下WFD ＜6 % ,TED ＜3 % ，整体误差在可接受范围内；n>25时WFD逐渐趋于一固定值，不再继续减小，这说明阶次增加已很难使FIR与IIR特性更加接近.综合考虑系统误差与实用性，本文认为“T-tube”模型特性可由15到25阶FIR函数近似，这将为MBSI算法的阶次估计提供依据.

图 6(Figure 6)

Fig. 6

图6 FIR与“T-tube”输出间的误差变化趋势 Fig. 6 The deviation trend of the FIR and‘T-tube’outputs a)—正常生理状态； (b)—血管收缩状态； (c)—血管舒张状态.

本文采用“黑箱”方法将各通路传递函数视为FIR函数，其阶次与系数均可由CR算法直接估计，大大降低了辨识难度. “T-tube”模型的验证结果表明，MBSI算法稳定性与准确性较好，能及时辨识心血管系统的不同状态.MBSI算法产生的收缩压误差小于4.6 mmHg,舒张压误差小于 4.2mmHg.

本文提出的MBSI算法无需建立先验函数或模型，实现了CAP的无创、个体化测量，但该算法应用于临床还有几个关键问题需要解决.论文采用的“黑箱”模型生理意义尚不明确，虽已证明脉搏波传导特性可由FIR函数描述，但其参数无法与心血管系统生理特性建立联系.盲辨识技术本身还有一些理论问题需要解决，如系统阶次确定方法、可辨识条件判定标准等.