刘文彦, 徐礼胜, 李宗鹏, 姜志豪
东北大学 中荷生物医学与信息工程学院, 辽宁 沈阳 110169
收稿日期：2018-09-11
基金项目：国家自然科学基金资助项目(61773110)；中央高校基本科研业务费专项资金资助项目(N161904002)。
作者简介：刘文彦(1991-)，女，山西大同人，东北大学博士研究生;
徐礼胜(1975-)，男，安徽安庆人，东北大学教授，博士生导师。

摘要：有创测量中心主动脉压力波形(central aortic pressure waveform, CAPW)存在一定的风险且花费较大.通常用外周动脉压力波形(peripheral arterial pressure waveform, PAPW)来代替CAPW, 而PAPW与CAPW有差异性.虽然盲系统辨识方法能够实时重建与动态辨识中心主动脉压力波形, 但该方法在模型定阶上还需改进, 因此本文提出了基于Hankel矩阵的秩定阶结合子空间方法用外周动脉压力波形重建中心主动脉压力波形.临床数据验证算法的结果表明，实测的与估计的CAPW的均方根误差约为5.87 mmHg, 波形匹配度约为74.70%;动物数据进一步验证算法, 结果表明实测的与估计的CAPW的均方根误差约为8.82 mmHg, 波形匹配度为54.57%.实验验证了该方法的有效性.
关键词：主动脉脉搏波外周动脉脉搏波盲系统辨识子空间法Hankel矩阵
Reconstruction of Aortic Pressure Wave Based on Subspace Algorithm
LIU Wen-yan, XU Li-sheng, LI Zong-peng, JIANG Zhi-hao
School of Sino-Dutch Biomedical & Information Engineering, Northeastern University, Shenyang 110169, China
Corresponding author: XU Li-sheng, E-mail: xuls@bmie.neu.edu.cn
Abstract: Invasive measurement of central aortic pressure waveform(CAPW)is risky and costly. Peripheral arterial pressure waveform(PAPW)has been used to substitute CAPW. However, there is a difference between PAPW and CAPW. Although the blind system identification technology can reconstruct in real-time and dynamically estimate CAPW, it needs to be improved in optimizing the model ordering. Accordingly, a rank-based ordering method was proposed based on Hankel matrix combining the subspace method to reconstruct the central aortic pressure waveform with two peripheral pressure waveforms. The results of the validation algorithm based on the clinical data showed that the root mean square error of the measured CAPW and the estimated CAPW is about 5.87 mmHg, and the waveform fit is about 74.70%. The results of the validation algorithm based on the collecting animal data showed that the root mean square error of the measured CAPW and the estimated CAPW is about 8.82 mmHg, and the waveform fit is about 54.57%. Thus the experiments verified the effectiveness of the method.
Key words: aortic pulse waveformperipheral arterial pulse waveformblind system identificationsubspace algorithmHankel matrix
随着社会快速发展, 城市的快节奏导致人们饮食习惯、生活方式等发生了巨大变化，加剧了心血管疾病的发病率及死亡率.心血管疾病是全球范围内导致人类死亡的主要原因, 今后10年心血管疾病患病人数仍处在快速上升趋势^[1].如何检测与预防心脑血管疾病的发病率及死亡率仍是全社会亟需解决的重大难题.中心主动脉压力波形(central aortic pressure waveform, CAPW)对评估心脑血管血流动力学有非常重要的价值.在主动脉根部测得的CAPW直接反映了心脑血管压力参数, 其可能成为心血管事件预测指标, 所以准确获得CAPW对检测与预防心脑血管疾病有重大意义^[2].目前中心主动脉压测量的“金标准”采用有创测量.由于有创测量创伤大、并发症多、费用昂贵等原因，该方法有一定的局限性, 近年来研究者们提出用无创方法来估计中心动脉压力^[3].无创测量主要是传递函数法和替代法, 传递函数利用外周动脉压力波形通过函数转换或模型间接推导CAPW, 该方法更具有实用性, 但传递函数法需要建立先验传递函数，而心血管系统个体间有差异性及时变性, 所以该方法在精度、实际应用上还需要进一步提高.替代法常用上臂肱动脉血压代替中心主动脉血压评估心血管系统功能特性^[4-6].颈动脉被认为接近中心主动脉, 常用颈动脉波形代替中心主动脉波形，可以很好地反映心脑血管的压力, 但颈动脉脉搏波无创采集需要借助超声测量血管壁直径的变化, 其操作复杂, 而且动脉系统中波的传播与反射现象导致不同部位压力波形的差异性, 故外周动脉压力波形(peripheral arterial pressure waveform, PAPW)不能直接代替CAPW^[7-8].近年来盲系统辨识技术在生物医学领域的广泛应用, Hahn等提出多路外周动脉脉搏波估计中心主动脉脉搏波的多通道盲系统辨识方法^[9].但目前盲系统辨识(blind system identification, BSI)算法在理论上与临床应用上还不够成熟, 仅仅根据已知外周动脉压力波形估计系统传递函数的阶数及参数等问题还有待解决.本文利用双通道盲系统辨识技术, 先用肱动脉与股动脉压力波形, 通过具体的盲系统辨识法获得信道矩阵，再构建Hankel矩阵来确定系统传递函数的阶数，并进一步辨识获得信道参数, 最后估计出中心主动脉压力波形.
1 数据与方法1.1 数据的采集1) 临床病人数据.由于实验条件限制及考虑病人手术时间等特殊情况, 同步采集三路有创数据比较困难.本文采集了25位病人的中心主动脉脉搏波和上肢肱动脉脉搏波压力波形数据.在患者知情情况下, 所有数据于中国医科大学附属第一医院心内科介入导管手术室采集, 采集的两路脉搏波数据是借助有创导管测量.受试者均是心脏介入手术的患者.患者的选择标准是不限制性别、年龄, 患有心血管疾病并需要进行手术的患者, 同时排除了休克等可能影响患者意识和血压测量的情况.
2) 动物实验数据.本文选择兔子作为实验动物, 选择标准是排除兔子本身心律失常、心室颤动等疾病, 体重为2.0~3.0 kg.本文的目的是通过有创导管法同步采集兔子的三路脉搏波, 具体测量部位是主动脉、股动脉、肱动脉.本文开始通过有创测量颈总动脉压力来代替主动脉压力, 在后续工作中进一步测量了主动脉压力.
1.2 多通道盲系统辨识算法模型盲系统辨识是指从系统已知的输出信号辨识出系统未知的信息, 主要包括系统传递函数的参数和系统输入信号.BSI广泛应用在不同的领域, 如数据通信、语音识别、地震勘察等，这些领域的共同难点是输入信号获得非常困难、代价大、不可测量.这引起了研究者关注并引出了大量的方法^[10-13], 如交叉相关(cross-relation, CR)法、两步极大似然(two-step maximum likelihood, TSML)法、子空间(subspace, SS)法, 这些方法都是基于二阶统计量的方法.根据研究比较^[14-15], SS是最优的一种方法.图 1是单输入多输出(single-input multiple-output, SIMO)系统模型.
图 1(Fig. 1)

图 1 SIMO系统模型Fig.1 Single-input multiple-output system model

SIMO系统表示如下:

(1)

式中:L为信道的阶次; M为信道的个数; n=0, 1, …, N-W, 其中N为输出数据的长度; W为窗参数, 它决定了每个输出向量y的长度.y定义如下:

(2)

(3)

输入序列s定义如式(4)所示:

(4)

为了推导方便, 假设噪声为加性高斯白噪声方差δ², 其他未知噪声同样满足子空间算法^[16], 下面为高斯白噪声向量:

(5)

(6)

构建信道传输矩阵H_M, M个子信道系数向量为

(7)

其中H_(i)表示如下:
式(1)改写为式(8):

(8)

1.3 子空间法子空间法是通过已知输出y(n)计算相关协方差矩阵R_y, 利用信号子空间与噪声子空间的正交关系求解.式(8)两边分别求相关协方差矩阵，得到式(9):

(9)

R_y的奇异值分解如式(10)所示:

(10)

令λ_i表示R_y的特征值, i=0, 1, …, MW－1，且有λ₀≥λ₁≥…≥λ_MW－1, 这里R_V=δ²I.传输矩阵H_M在盲系统辨识中起到关键作用.只有当该矩阵为列满秩才能进行盲系统辨识, 因此H_M列满秩也作为子空间法的假设条件之一, 其实在实际应用中也满足此条件^[16].由上面关系可得:λ_i>δ², i=0, …, W+L－1;λ_i=δ², i=W+L, …, MW－1.U_S=[u₀, …, u_W+L－1]是λ₀, …, λ_W+L－1对应的特征向量组成的向量空间，称为信号子空间, 记为Span{U_S}.U_V=[u_W+L, …, u_MW]是λ_W+L, …, λ_MW对应的特征向量组成的向量空间，称为信号噪声子空间, 记为Span{U_V}.而H_M的列向量展开同样是信号子空间, Span{U_S}=Span{H_M}.信号子空间与噪声子空间的正交关系是进行盲系统辨识的重要基础，可得到式(11):

(11)

令g_i(0≤i≤MW－W－L－1)表示U_V的列向量, 其可分割成M个小的向量, 即g_i=[g_i1^T, …, g_iM^T]^T, G_ij是利用g_ij形成Toeplitz结构的矩阵, 其中第一行为[g_ij(0), g_ij(1), …, g_ij(W－1), ], 第一列为, G_i=[G_i1^T, …, G_iM^T]^T∈C^{M(L+1)×(W+L)}.考虑到输出系统中会有噪声, 因此这里采用线性最小均方误差估计准则来求解式(11).记损失函数Q为

(12)

SS法就是求q在的约束条件下, 对q特征分解, 求最小特征值对应的特征向量求出.考虑到这里求出的是一个估计量, 对实际的h有一个尺度变化.本文对h进行了约束, 令其2范数等于1, 并估计出2个通道参数及中心主动脉压力波形.
1.4 模型的阶数估计对于盲系统辨识过程来说, 只有输出数据已知情况下, 利用盲系统辨识算法可获得有限冲激响应序列, 为了进一步把有限冲激响应序列转变成参数模型, 需要确定模型的阶数.本文利用有限冲激响应构建Hankel矩阵, 求矩阵的秩来判断模型的阶数^[17].本文在Hankel矩阵行列式的平均比值法的基础上, 为了进一步提高精度用Hankel矩阵的行列式绝对值比值法来判断系统阶次.假设盲系统辨识中有限冲激响应序列记作h(1), h(2), …, h(L)，构造如下Hankel矩阵:

(13)

其中:l决定Hankel矩阵的维数; k在1到(L-2l+2)之间任意选择, 它决定用哪些冲激响应序列组成Hankel矩阵.
1) FIR序列不含噪声情况.假设h(1), h(2), …, h(L)不含噪声, 则判断盲系统辨识系统模型的阶数估计如下:按式(13)构造的Hankel矩阵H(l, k), 对给定的l, 计算k取1到L-2l+2时Hankel矩阵的行列式的值.如果l从1逐一增加到n₀, 对所有k, 都有det[H(l, k)]≠0;而l从l增加到n₀+1后, 对所有的k, det[H(l, k)]=0, 说明Hankel矩阵在l=n₀+1处, 由非奇异矩阵变成奇异矩阵, 则可以判断模型的阶数为n₀.
2) FIR序列含噪声情况.假设h(1), h(2), …, h(L)含噪声, 因此l从1增加到n₀+1后, 对所有的k, Hankel矩阵的行列式都不会为零, 因此不能用无噪声情况去判断模型的阶数.
对模型阶数确定方法的改进, 用Hankel矩阵行列式的平均值比值D_l来观察Hankel矩阵是否由非奇异矩阵变为奇异矩阵.D_l定义如下:

(14)

当l从1开始逐一增加时, 不断计算D_l的值，可取D_l到达最大值的l作为模型的阶数.
考虑到式(13)中, Hankel矩阵行列式有正有负, 计算行列式平均值时出现正负抵消情况, 为了提高辨识精度进一步改进算法, 用Hankel矩阵行列式绝对值的平均值比值变化来确定模型阶数.

(15)

当l从1开始逐一增加时, 不断计算E_l的值，可取E_l到达最大值的l作为模型的阶数.
1.5 对方法准确度评价为了定量评价盲系统辨识法的性能, 需要一些合适的性能评价指标多角度判断.常用的评价准则有基于信号和系统矩阵的评价准则^[18].本文利用基于信号本身的一些特性作为评价的准则.
1) 波形匹配度.波形匹配度(waveform fit, WF)大小反映原信号与盲系统辨识法估计的信号波形之间的相关系数的大小.WF越趋近于1, 表明盲系统辨识法估计的信号波形与实测的主动脉波形相似性越高、辨识效果越好.

(16)

式中：y(i)是盲系统辨识法估计出的中心主动脉信号；s(i)是实测的中心主动脉信号；N是采样点数.
2) 均方根误差.均方根误差(root mean square error, RMSE)用来衡量实测中心主动脉波形与盲系统辨识估计中心主动脉波形之间的偏差，

(17)

2 结果通过“T-tube Load”模型, 利用实测的中心主动脉压力波形数据作为模型的输入估计出两路外周脉搏波压力波形.如图 2所示, 仿真得到的两路外周压力脉搏波波形明显比实测的中心主动脉压力脉搏波波形的收缩压及平均压要大, 符合动脉脉搏的生理性放大效果.通过仿真模型可以得到两路外周脉搏波压力波形数据来验证盲系统辨识法.
图 2(Fig. 2)

图 2 实测的中心主动脉压力波形与“T-tube Load”模型估计的肱动脉、股动脉压力波形Fig.2 Measured central aortic pressure waveform and estimated pressure waveform of brachial artery and femoral artery using the "T-tube Load" model

本文比较了三种盲系统辨识算法, 在模型定阶上都用Hankel矩阵的秩来确定.如图 3所示, 估计的CAPW与实测的CAPW之间的关系表明SS法最好、其次TSML法、最后为CR法.本文的定阶方法对于CR法、TSML法可能还存在问题.
图 3(Fig. 3)

图 3 三种算法估计与实测的CAPW之间的关系Fig.3 Relationship of estimated and measured CAPW using three algorithms (a)—RMSE; (b)—WF.

本文最终选择用SS法来估计中心主动脉波形, 其中系统的阶次估计用Hankel矩阵的秩来定阶, 通过计算Hankel矩阵行列式绝对值的变化来确定系统传递函数的阶次.
本文通过Hankel矩阵的秩定阶结合子空间法, 如图 4所示, 利用仿真得到的两路脉搏波压力波形来估计出中心主动脉压力波形与实测的中心主动脉压力波形的RMSE约为5.87 mmHg, WF约为74.70%.如图 5所示, 通过动物实验实测的主动脉压力波形与估计的中心主动脉压力波形之间的RMSE平均值为8.82 mmHg, WF平均值为54.57%.
图 4(Fig. 4)

图 4 SS法估计与实测的CAPWFig.4 Estimated and measured CAPW using the SS algorithm

图 5(Fig. 5)

图 5 SS法估计与实测的CAPWFig.5 Estimated and measured CAPW using the SS algorithm

3 结语通过“T-tube Load”模型, 利用有创中心主动脉压力波形来估计两路外周脉搏波压力波形验证子空间法.通过仿真数据验证SS法可行性, 然后进一步采集动物数据来验证SS法的可行性.本文在模型定阶上解决了盲系统辨识算法对阶数的确定问题, 但实测的与估计的主动脉压力波形精度还有待进一步提高, 将来的工作需要继续改进算法.本文中提到的双通道盲系统辨识算法有几个关键问题, 选择哪两路外周脉搏波组合来估计中心主动脉脉搏波是首要任务, 在将来研究中要考虑外周脉搏波的选择; 盲系统辨识算法中对模型阶数辨识方法, 还需要不断的改进.
参考文献

[1]	陈伟伟, 高润霖, 刘力生, 等. 《中国心血管病报告2017》概要[J].中国循环杂志, 2018, 33(1): 3–6. ( Chen Wei-wei, Gao Run-lin, Liu Li-sheng, et al. China cardiovascular disease report 2017[J].Chinese Circulation Journal, 2018, 33(1): 3–6.)
[2]	Suleman R, Padwal R, Hamilton P, et al. Association between central blood pressure, arterial stiffness, and mild cognitive impairment[J].Clinical Hypertension,, 2017, 23(1): 2–6.
[3]	Hajjar I, Goldstein F C, Martin G S, et al. Roles of arterial stiffness and blood pressure in hypertension-associated cognitive decline in healthy adults[J].Hypertension, 2015, 67(1): 171–175.
[4]	Cai T Y, Qasem A, Ayer J G, et al. Central blood pressure in children and adolescents:non-invasive development and testing of novel transfer functions[J].Journal of Human Hypertension, 2017, 31(12): 831–837.DOI:10.1038/jhh.2017.59
[5]	Gao M, Rose W C, Fetics B, et al. A simple adaptive transfer function for deriving the central blood pressure waveform from a radial blood pressure waveform[J].Scientific Reports, 2016, 6(1): 1–9.DOI:10.1038/s41598-016-0001-8
[6]	Daisuke S, Eiichiro Y, Tomoko T, et al. The accuracy of central blood pressure waveform by novel mathematical transformation of non-invasive measurement[J].International Journal of Cardiology, 2015, 189(1): 244–246.
[7]	Lee J C, Ghasemi Z, Kim C S, et al. Investigation of viscoelasticity in the relationship between carotid artery blood pressure and distal pulse volume waveforms[J].IEEE Journal of Biomedical & Health Informatics, 2017, 22(2): 460–470.
[8]	Kips J G, Schutte A E, Vermeersch S J, et al. Comparison of central pressure estimates obtained from SphygmoCor, Omron HEM-9000AI and carotid applanation tonometry[J].Journal of Hypertension, 2011, 29(6): 1115–1120.DOI:10.1097/HJH.0b013e328346a3bc
[9]	Hahn J O, Reisner A T, Asada H H. Blind identification of two-channel ⅡR systems with application to central cardiovascular monitoring[J].Journal of Dynamic Systems Measurement & Control, 2009, 131(5): 545–553.
[10]	Qiu W, Hua Y.Performance comparison of three methods for blind channel identification[C]// IEEE International Conference on Acoustics.[S.l.]: ICASSP, 1996: 2423-2426.
[11]	Xu G, Liu H, Tong L, et al. A least-squares approach to blind channel identification[J].IEEE Transactions on Signal Processing,, 1995, 43(12): 2982–2993.DOI:10.1109/78.476442
[12]	Hua Y B.Fast maximum likelihood for blind identification of multiple FIR channels[C]//Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers.Pacific Grove, 1994: 471-478.
[13]	Patel A, Li J, Finegan B, et al. Aortic pressure estimation using blind identification approach on single input multiple output non-linear Wiener systems[J].IEEE Transactions on Biomedical Engineering, 2018, 65(6): 1193–1200.DOI:10.1109/TBME.2017.2688425
[14]	Mayyala Q, Meraim K A, Zerguine A. Structure-based subspace method for multi-channel blind system identification[J].IEEE Signal Processing Letters, 2017, 24(8): 1183–1187.DOI:10.1109/LSP.2017.2715418
[15]	Moulines E, Duhamel P, Cardoso J F, et al. Subspace methods for the blind identification of multichannel FIR filters[J].IEEE Transactions on Signal Processing, 1995, 43(2): 516–525.DOI:10.1109/78.348133
[16]	Diamantaras K I, Papadimitriou T. An efficient subspace method for the blind identification of multichannel FIR systems[J].IEEE Transactions on Signal Processing, 2008, 56(12): 5833–5839.DOI:10.1109/TSP.2008.929132
[17]	方崇智. 过程辨识[M]. 北京: 清华大学出版社, 1998. ( Fang Chong-zhi. Process identification[M]. Beijing: Tsinghua University Press, 1998.)
[18]	马春华.基于子空间的盲辨识和盲均衡算法[D].西安: 西安电子科技大学, 2005. ( Ma Chun-hua.Subspace-based methods for blind identification and blind equalization[D].Xi'an: Xidian University of Electronic Technology, 2005.http://cdmd.cnki.com.cn/article/cdmd-10701-2005043169.htm)