心阻抗血流图（ICG）在心血管疾病患者心功能评价中发挥重要作用。针对ICG信号的测量容易受运动伪迹的干扰问题，本文提出了两步频域集合经验模态分解（EEMD）和典型相关性分析（CCA）的去噪方法。首先，将ICG信号、心电图（ECG）分别与运动信号做第一步频域EEMD-CCA，将相关性系数高的成分置零来抑制主要运动伪迹。其次，将得到的ECG和ICG信号做第二步频域EEMD-CCA，获取这两组生理信号之间的共同成分来进一步去噪。最后，利用这些共同成分来重构ICG信号。本研究招募了30名志愿者参与测试，结果表明，使用本文提出的去噪方法后，ICG信号的质量有大幅度的提升，可为后续的心血管疾病诊断和分析提供支撑。
引用本文： 解尧, 杨东, 余洪龙, 解启莲. 基于两步频域集合经验模态分解和典型相关性分析的运动心阻抗血流图信号的去噪方法研究. 生物医学工程学杂志, 2024, 41(5): 986-994. doi: 10.7507/1001-5515.202210059 复制
0 引言
每搏量（stroke volume，SV）是指一个心动周期内从左心室射出的血液量。心排量（cardiac output，CO）是由SV与心率（heart rate，HR）相乘得到的，代表每分钟从心脏排出的血液量[1]。无创、连续地估算SV和CO，对于评估心衰、冠心病以及其他心血管疾病意义重大[2]。测量SV的经典方法是指示剂稀释法，主要包括弗克（Fick）法、染料稀释法和热稀释法[3-5]。然而，由于具有有创和非连续性的测量特点，使得指示剂稀释法的使用场景受到严重限制。食管超声心动图和多普勒超声心动图可以无创测量SV[6]，但仍然是以非连续性的方式测量，且测量结果的准确性也依赖于操作者的经验。心阻抗血流图（impedance cardiography，ICG）是一种无创、连续测量SV和CO的方法，它通过监测胸部阻抗，并将阻抗周期性变化与心脏周期性搏动联系起来[7]。SV和CO的测量依赖于算法对ICG信号波形中特征点的准确识别。如图1所示，关键特征点包括Q、B、C和X，分别代表心室除极、主动脉瓣的打开时刻、最大射血速度时刻和主动脉瓣的关闭时刻[8]。然而，当身体运动和呼吸产生运动伪迹时，这些特征点就难以被准确检测。
图1
ICG波形的特征位置关系
Figure1.
Location relationship of ICG features
现有去除ICG信号运动伪迹的方法有多种，各有特点。其中，带通滤波，可以去除频带外的噪声成分[9]，但在处理与目标信号频谱重叠的伪迹时具有一定的局限性。自适应滤波，是另一种常用的去噪方法，可以自动更新滤波器的权重，以适应不同强度的噪声[10-11]。虽然相关研究已经将自适应滤波技术应用于ICG信号处理，但通常很难确定各种伪迹来源的滤波器权重。小波变换，是数字信号的时间尺度表示，已广泛应用于ICG信号去噪；基于小波变换的方法，其难点是选择最优母小波[12-13]，它常常对去噪效果起着决定性作用。集成平均，是通过对每个信号的数字化样本进行平均，从而减少信号单拍波动的影响[14-15]；尽管这种方法可以消除由情绪紧张等引起的单拍伪迹，但在运动强度较大的情况下，去噪效果不佳。
基于上述研究的成果与不足，本文提出基于两步频域集合经验模态分解（ensemble empirical mode decomposition，EEMD）和典型相关性分析（canonical correlation analysis，CCA）的去噪方法[16-17]。EEMD可将单通道数据分解为多通道数据；CCA可用于计算两组变量之间的整体相关性。由于CCA需要多通道参数输入，通常需要EEMD配合使用，即EEMD-CCA，故而本文提出的方法可简称为两步频域EEMD-CCA的去噪方法。第一步频域EEMD-CCA是考虑到运动伪迹来源于运动[18-19]，那么ICG信号与运动信号之间，以及心电图（electrocardiogram，ECG）信号与运动信号之间，最相关的部分应为运动伪迹；第二步频域EEMD-CCA是分析了ICG信号和ECG信号的来源，由于两者都反映心脏周期性舒张和收缩，它们之间最相关部分应为ICG信号的主要成分[20]。
1 实验材料与设备
本文验证数据集来源于30名受试者（17名男性，13名女性，年龄在18～55岁之间，身高在155～183 cm），分别在静息状态和运动状态下测量采集得到。试验获得安徽医科大学附属第四医院伦理委员会批准（伦理批号：PJ2021-006-01），同时所有受试者都签署知情同意书。试验的数据处理分析已经获得受试者同意。受试者入组标准：① 年龄18周岁以上，性别不限；② 受试者需要进行CO检测或者自愿参加CO检测；③ 受试者能够理解研究目的，并签署知情同意书。排除标准：① 胸部创伤、胸部手术患者；② 电极片皮肤过敏的患者；③ 放置心脏起搏器的患者；④ 体外循环辅助患者等。
所有信号采集持续时间为150 s，如图2所示，受试者坐在扶手椅上（采集30 s）作为静息状态，骑乘动感单车（采集120 s）作为运动状态。ICG信号采集电极，ECG信号采集电极和加速度（acceleration，ACC）测量模块（简称：加速度计）的位置分别用绿色、红色和橙色框标记。通过无创CO监测仪（TL-NiCON-200，安徽通灵仿生科技有限公司，中国），同时获得了ICG信号、ECG信号和加速度信号。ICG信号是根据Kubicek等[7]描述的四电极阵列获取。将电流电极放置在颈部底部，接地电极靠近剑胸交界下方处，两个电压电极分别与电流、接地电极相邻粘贴；注入的恒流源为2 mA和50 kHz。ECG信号是通过II导联阵列粘贴的两个电极获得，ECG信号接地电极与 ICG信号共用。无创CO监测仪集成了加速度传感器，可采集x、y、z三个方向的实时加速度信号。算法的进一步分析和研究在商业数学软件MATLAB R2016a（MathWorks Inc.，美国）上完成。
图2
受试者的ICG和ECG信号的采集场景
Figure2.
Acquisition scenario of ICG and ECG signals of subjects
2 去噪方法研究
2.1 方法介绍
本文算法流程主要包括4大部分：预处理、第一步频域EEMD-CCA、第二步频域EEMD-CCA、ICG信号重构，如图3所示。预处理，对生理信号进行带通滤波，融合运动信号。第一步频域EEMD-CCA，旨在抑制ICG信号和ECG信号中的主要运动伪迹。第二步频域EEMD-CCA，是为了获取ICG信号和ECG信号主要成分来进一步去除噪声。最后，重构ICG信号，通过两个信号最相关部分来完成重构。
图3
ICG信号的去噪流程
Figure3.
Denoising process of ICG signal
2.1.1 预处理
首先将ECG信号、ICG信号和加速度信号切片成不重叠的片段，每个片段切片长度为10 s。ICG信号的主要频段在 [0.5, 4.0] Hz的频带内[21]，考虑到运动过程中心率的增加，ICG信号的主成分频谱将变宽。因此，ICG 信号和加速度信号的截止频率设置为 [0.7, 6.0] Hz，ECG 信号的截止频率设置为 [0.7, 25.0] Hz。对x、y、z方向的加速度信号做快速傅里叶变换，将变换后获得的幅度部分分别标记为向量Fx(k)、Fy(k)和Fz(k)。选择三个运动方向上频谱幅度最大值作为合并的运动加速度信号能有效地表征运动伪迹[22]，如式(1)所示：
 |
其中，FM(k)为合并的加速度信号的频谱向量，即运动信号的频谱向量。
2.1.2 两步频域
根据CCA执行过程中的要求[23]，需要将单通道ECG信号、ICG信号和加速度信号分解为多通道信号集，在这里使用EEMD将单通道的ICG信号、ECG信号以及加速度信号分解成多通道的本征模态函数（intrinsic mode function，IMF）集合。EEMD通过将真实的IMF定义为试验集合的平均值，并且每个信号由原始信号加上一个有限幅度的白噪声组成，从而消除了模态混叠的困境[24-25]。
EEMD将信号x(t)（ICG信号、ECG信号或加速度信号）与白噪声n(t)相加，生成新的信号y(t)，再进一步分解为IMF，具体如式(2)和式(3)所示：
 |
 |
其中，zi(t)为x(t)信号分解的第i个IMF，K为得到的IMF总数，α为白噪声的幅值，rn(t)为EEMD分解后的残余信号。
由于ICG信号主频率范围为［0.7, 6.0］ Hz，本研究只选择主频率在这个范围内的IMF。为了保证后续CCA处理的性能，要求选择所选IMF的通道数至少为5个，如果实际符合要求的IMF少于5个，保留主频率邻近的IMF作为补充。
将所有被选中的IMF做快速傅里叶变换，并取其幅度部分分别记为：s1, ICG, ···, sN, ICG、 s1, ECG, ···, sN, ECG、s1, ACC, ···, sN, ACC，其中N为被选中的IMF个数，N ≤ K。定义ICG、ECG以及加速度信号的特征矩阵，SICG = ［s1, ICG, ···, sN, ICG］，SECG = ［s1, ECG, ···, sN, ECG］，SACC =［s1, ACC, ···, sN, ACC］。SICG、SECG以及SACC是频域CCA处理的输入信号，两步频域CCA主要是为了计算SICG、SECG以及SACC之间相关性，再根据相关性系数进一步处理带有运动伪迹的ICG信号。
假定CCA的典型相关变量
和
可以用目标信号Sa和Sb来表示：
、 
，其中，w1和w2为权值向量，μa和νb为典型相关变量[23]，最终CCA的目标函数定义如式(4)所示：
 |
其中，r为典型相关系数向量，
和
分别为Sa和Sb的自协方差矩阵，
为Sa和Sb的交叉协方差矩阵。假设得到N个典型相关系数，用向量r来代表典型相关系数向量，r = ［r1, r2, ···, rk, ···, rN］；其中，0 ≤ ri ≤ 1。A和B的N个典型相关变量分别记录在矩阵Ua和Vb，相应的权值向量分别存储在矩阵Ca和Db中。
对于第一步频域CCA，计算ICG和加速度信号分解的IMF之间的频域相关性系数，以及ECG和加速度信号分解的IMF之间的频域相关性系数。这里认为相关性系数较大的部分是运动引起的噪声，因此将相关性系数高（ri > 0.8）的成分对应的典型相关变量置0，来抑制运动伪迹。
对于第二步频域CCA，计算ICG和ECG信号分解的IMF之间的频域相关性系数。此时，认为相关性系数大的部分是两个信号的主频成分，因此将相关性低（ri < 0.9）的成分对应的典型相关变量置0，来进一步去除运动伪迹。
2.1.3 重构ICG信号
重构ICG信号， ICG频谱可以通过保留的权值系数矩阵CICG和对应的部分分量被置为0的典型相关矩阵UICG进行重构，如式(5)所示：
 |
其中，HICG = ［h1, ICG, ···, hN, ICG］是ICG信号重构的频谱，hi, ICG是ICG信号第i个IMF去噪后的频谱向量。
因此，ICG信号的第i个IMF去噪后的时域信号zi, ICG就可以利用频谱向量hi, ICG信息和对应的相位信息，再通过傅里叶逆变换计算得出，如式(6)所示：
 |
其中，zi, ICG为第i个IMF去噪后的时域信号，M为傅里叶逆变换时的点数，π为圆周率，k为频域信号的频率索引，j为虚数单位，n为时域信号的时间索引。最后，将所有得到的时域IMF相加获得最终重构的ICG信号，如式（7）所示：
 |
其中，yICG为最终重构的ICG信号，N为IMF的个数，n为时域信号的时间索引。
2.2 信号质量的检验指标
为了评估ICG信号去噪的有效性，需要对信号质量做进一步的评估，这是因为ICG的信号质量影响特征点检测的准确性，进而影响相应生理参数的计算准确性。由于缺乏评估ICG信号的量化指标，本文给出了通过ICG信号得到的相关生理参数的合理范围，如表1所示，并通过该范围限制来判定信号在去噪后的效果。一个有效的心拍应满足以下条件：① 射血前期（pre-ejection period，PEP）在50～160 ms的范围内；② 左室射血时间（left ventricular ejection time，LVET）在160～380 ms的范围内；③ 阻抗微分最大值在0.4～3.0 Ω/s的范围内；④ HR在30～180 beat/min的范围内[26]。
表1
ICG信号特征参数的正常值范围
Table1.
Normal range of ICG signal characteristic parameters
为了进一步验证信号质量，本文增加了更多细化的评价指标，主要包括：① 有效节拍（valid beat，VB）；② 节拍贡献因子（beat contribution factor，BCF），其计算方式为VB数目与总节拍数的比值，BCF越大说明信号质量越好；③ CC，为ICG信号中两个相邻C点的时间间隔，反映的是心动周期；CC计算出的心动周期与实际心动周期差值越小，信号质量越高；④ PEP变异性和LVET变异性，即在连续ICG数据的片段内（10 s），确认每个片段内测量的PEP和LVET的差异是否较小，亦即方差是否较小。
这些评价指标的测量依赖对ICG信号C点、B点以及X点等特征点的识别。ICG的C点是最大射血速度的时刻，是ICG在一个心动周期的幅度值最大的点[27-28]。ICG的B点是主动脉瓣的打开时刻，是幅值为C点幅值的0.15倍的特征点[29-31]。ICG的X点是主动脉瓣的关闭时刻，是ICG在一个心动周期的幅度值最小的点[27]。ICG信号的特征点检测方法可以总结如下：
（1）搜索ECG信号，识别特征点Q、R、T。
（2）扫描对应时刻的ICG信号，从ECG的R点对应时刻开始的三分之一个心动周期内，搜索ICG信号幅度的最大值，这个最大值对应特征点记为C点。
（3）扫描对应时刻的ICG信号，从ECG的R点对应时刻开始到ICG的C点对应时刻，找到幅值为C点幅值的0.15倍的特征点，并将其记为B点。
（4）扫描对应时刻的ICG信号，从ECG的T点对应时刻开始的三分之一心动周期内，将幅度最小值的特征点标记为X点。
（5）重复上述步骤，直到ICG信号结束。
3 结果与分析
3.1 结果
如图4所示，本文展示了ICG信号去噪效果和特征点识别检测，同时对ICG信号进行归一化，以便清晰地展示原始ICG和去噪后的ICG之间的差异。ICG信号的不同特征点用不同的颜色和标记来表示，可以清晰发现去噪后的信号获得了更精确的特征点。如表2所示，为30名受试者的测量结果。其中，静息测量场景共1 068次心拍；去噪后，ICG的VB数量从976次增加到1 055次，ICG的BCF也从原来的91.4%上升到98.7%。在运动测量场景中也有类似的改善，特别是ICG的BCF由78.4%提高到97.5%。总体而言，去噪前ICG信号的平均BCF约为80.1%，去噪后ICG的平均BCF约为97.7%。
图4
ICG信号去噪效果以及特征点识别
Figure4.
Denoising effect of ICG signal and feature point detection
表2
静止和运动两种状态下去噪效果的结果
Table2.
Results of the denoising effect in rest and motion state
由于ECG信号的R波特征明显，即使在运动伪迹干扰下也较为容易识别和测量，可以作为ICG信号的一个参考。如图5所示，去噪前估算出的心动周期与实际心动周期之间存在较大的误差值，尤其是在运动状态下。相比之下，去噪后估算出的心动周期和实际心动周期之间误差值明显减小，同时估计值与真实值的差值落在95%一致性区间内的数量更多。综上所述，在整个研究场景中（包括静息状态和运动状态），与去噪前相比，去噪后估算出的心动周期与实际心动周期更加一致，间接说明了ICG信号质量在去噪后得到提高。
图5
不同状态下心动周期的真实值与估算值的一致性
Figure5.
The consistence between the true value and the estimated value of the cardiac cycle in different states
如图6所示，对去噪后的ICG信号切片分组后，每个切片内数据获取的多个LVET值的方差均值明显低于去噪前。对于PEP值也有类似的结果，去噪后每个信号切片测量的PEP值的方差明显降低，间接说明去噪后ICG信号质量的提高。
图6
ICG信号在去噪前后，PEP和LVET参数的标准差的变化
Figure6.
Changes in the standard deviation of PEP and LVET parameters detected before and after ICG signal denoising
3.2 消融实验
上述实验结果验证了两步频域EEMD-CCA去噪方法对带有运动伪迹的ICG信号去噪的有效性。由于提出的去噪方法需要两步频域EEMD-CCA，本文通过消融实验来验证方法中每一步的必要性。第一步频域EEMD-CCA消除了噪声引起的主要频率成分，而第二步频域EEMD-CCA获取生理信号频谱中最相关的部分，进一步去除了运动产生的伪迹。如图7所示，单独使用第一步频域EEMD-CCA的去噪效果低于使用两步频域EEMD-CCA。如图8所示，单独使用第二步频域EEMD-CCA去噪效果也低于使用两步频域EEMD-CCA。因此，两步频域EEMD-CCA比仅仅使用其中一步频域EEMD-CCA在提高ICG质量方面更有效，验证了去噪方法的有效性。
图7
单独使用第一步频域EEMD-CCA与使用两步频域EEMD-CCA方法的去噪结果对比
Figure7.
Comparison of denoising results by using only the first step EEMD-CCA and the two-step EEMD-CCA method
图8
单独使用第二步频域EEMD-CCA方法与使用两步频域EEMD-CCA方法的去噪结果对比
Figure8.
Comparison of denoising results by using only the second step EEMD-CCA and the two-step EEMD-CCA method
4 总结
由于ICG信号的有效去噪有利于进一步的生理参数检测和心血管疾病诊断，故本文介绍了一种两步频域EEMD-CCA的方法对ICG信号进行去噪。本文方法在静息和运动两种状态下对30名受试者进行了评估。研究结果表明，去除运动伪迹后，ICG信号质量明显提高。
本文研究验证了所提去噪方法在静息和运动状态场景下的有效性，后续还需要进一步检验去噪方法在其他运动场景中的表现，如跑步、打球等运动幅度更大的场景下的有效性；同时，对ICG信号去噪只是心功能监测的一部分，后续还需要结合胸阻抗模型进一步验证相关心功能参数的准确性，如是否有利于对SV、CO计算准确度的提高。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：解尧负责软硬件平台搭建以及论文撰写、修订；杨东、余洪龙负责数据记录与分析；解启莲负责研究计划安排以及论文审阅修订。
伦理声明：本研究通过了安徽医科大学附属第四医院伦理委员会的审批（批号：PJ2021-006-01）。
0 引言
每搏量（stroke volume，SV）是指一个心动周期内从左心室射出的血液量。心排量（cardiac output，CO）是由SV与心率（heart rate，HR）相乘得到的，代表每分钟从心脏排出的血液量[1]。无创、连续地估算SV和CO，对于评估心衰、冠心病以及其他心血管疾病意义重大[2]。测量SV的经典方法是指示剂稀释法，主要包括弗克（Fick）法、染料稀释法和热稀释法[3-5]。然而，由于具有有创和非连续性的测量特点，使得指示剂稀释法的使用场景受到严重限制。食管超声心动图和多普勒超声心动图可以无创测量SV[6]，但仍然是以非连续性的方式测量，且测量结果的准确性也依赖于操作者的经验。心阻抗血流图（impedance cardiography，ICG）是一种无创、连续测量SV和CO的方法，它通过监测胸部阻抗，并将阻抗周期性变化与心脏周期性搏动联系起来[7]。SV和CO的测量依赖于算法对ICG信号波形中特征点的准确识别。如图1所示，关键特征点包括Q、B、C和X，分别代表心室除极、主动脉瓣的打开时刻、最大射血速度时刻和主动脉瓣的关闭时刻[8]。然而，当身体运动和呼吸产生运动伪迹时，这些特征点就难以被准确检测。
图1
ICG波形的特征位置关系
Figure1.
Location relationship of ICG features
现有去除ICG信号运动伪迹的方法有多种，各有特点。其中，带通滤波，可以去除频带外的噪声成分[9]，但在处理与目标信号频谱重叠的伪迹时具有一定的局限性。自适应滤波，是另一种常用的去噪方法，可以自动更新滤波器的权重，以适应不同强度的噪声[10-11]。虽然相关研究已经将自适应滤波技术应用于ICG信号处理，但通常很难确定各种伪迹来源的滤波器权重。小波变换，是数字信号的时间尺度表示，已广泛应用于ICG信号去噪；基于小波变换的方法，其难点是选择最优母小波[12-13]，它常常对去噪效果起着决定性作用。集成平均，是通过对每个信号的数字化样本进行平均，从而减少信号单拍波动的影响[14-15]；尽管这种方法可以消除由情绪紧张等引起的单拍伪迹，但在运动强度较大的情况下，去噪效果不佳。
基于上述研究的成果与不足，本文提出基于两步频域集合经验模态分解（ensemble empirical mode decomposition，EEMD）和典型相关性分析（canonical correlation analysis，CCA）的去噪方法[16-17]。EEMD可将单通道数据分解为多通道数据；CCA可用于计算两组变量之间的整体相关性。由于CCA需要多通道参数输入，通常需要EEMD配合使用，即EEMD-CCA，故而本文提出的方法可简称为两步频域EEMD-CCA的去噪方法。第一步频域EEMD-CCA是考虑到运动伪迹来源于运动[18-19]，那么ICG信号与运动信号之间，以及心电图（electrocardiogram，ECG）信号与运动信号之间，最相关的部分应为运动伪迹；第二步频域EEMD-CCA是分析了ICG信号和ECG信号的来源，由于两者都反映心脏周期性舒张和收缩，它们之间最相关部分应为ICG信号的主要成分[20]。
1 实验材料与设备
本文验证数据集来源于30名受试者（17名男性，13名女性，年龄在18～55岁之间，身高在155～183 cm），分别在静息状态和运动状态下测量采集得到。试验获得安徽医科大学附属第四医院伦理委员会批准（伦理批号：PJ2021-006-01），同时所有受试者都签署知情同意书。试验的数据处理分析已经获得受试者同意。受试者入组标准：① 年龄18周岁以上，性别不限；② 受试者需要进行CO检测或者自愿参加CO检测；③ 受试者能够理解研究目的，并签署知情同意书。排除标准：① 胸部创伤、胸部手术患者；② 电极片皮肤过敏的患者；③ 放置心脏起搏器的患者；④ 体外循环辅助患者等。
所有信号采集持续时间为150 s，如图2所示，受试者坐在扶手椅上（采集30 s）作为静息状态，骑乘动感单车（采集120 s）作为运动状态。ICG信号采集电极，ECG信号采集电极和加速度（acceleration，ACC）测量模块（简称：加速度计）的位置分别用绿色、红色和橙色框标记。通过无创CO监测仪（TL-NiCON-200，安徽通灵仿生科技有限公司，中国），同时获得了ICG信号、ECG信号和加速度信号。ICG信号是根据Kubicek等[7]描述的四电极阵列获取。将电流电极放置在颈部底部，接地电极靠近剑胸交界下方处，两个电压电极分别与电流、接地电极相邻粘贴；注入的恒流源为2 mA和50 kHz。ECG信号是通过II导联阵列粘贴的两个电极获得，ECG信号接地电极与 ICG信号共用。无创CO监测仪集成了加速度传感器，可采集x、y、z三个方向的实时加速度信号。算法的进一步分析和研究在商业数学软件MATLAB R2016a（MathWorks Inc.，美国）上完成。
图2
受试者的ICG和ECG信号的采集场景
Figure2.
Acquisition scenario of ICG and ECG signals of subjects
2 去噪方法研究
2.1 方法介绍
本文算法流程主要包括4大部分：预处理、第一步频域EEMD-CCA、第二步频域EEMD-CCA、ICG信号重构，如图3所示。预处理，对生理信号进行带通滤波，融合运动信号。第一步频域EEMD-CCA，旨在抑制ICG信号和ECG信号中的主要运动伪迹。第二步频域EEMD-CCA，是为了获取ICG信号和ECG信号主要成分来进一步去除噪声。最后，重构ICG信号，通过两个信号最相关部分来完成重构。
图3
ICG信号的去噪流程
Figure3.
Denoising process of ICG signal
2.1.1 预处理
首先将ECG信号、ICG信号和加速度信号切片成不重叠的片段，每个片段切片长度为10 s。ICG信号的主要频段在 [0.5, 4.0] Hz的频带内[21]，考虑到运动过程中心率的增加，ICG信号的主成分频谱将变宽。因此，ICG 信号和加速度信号的截止频率设置为 [0.7, 6.0] Hz，ECG 信号的截止频率设置为 [0.7, 25.0] Hz。对x、y、z方向的加速度信号做快速傅里叶变换，将变换后获得的幅度部分分别标记为向量Fx(k)、Fy(k)和Fz(k)。选择三个运动方向上频谱幅度最大值作为合并的运动加速度信号能有效地表征运动伪迹[22]，如式(1)所示：
|
其中，FM(k)为合并的加速度信号的频谱向量，即运动信号的频谱向量。
2.1.2 两步频域
根据CCA执行过程中的要求[23]，需要将单通道ECG信号、ICG信号和加速度信号分解为多通道信号集，在这里使用EEMD将单通道的ICG信号、ECG信号以及加速度信号分解成多通道的本征模态函数（intrinsic mode function，IMF）集合。EEMD通过将真实的IMF定义为试验集合的平均值，并且每个信号由原始信号加上一个有限幅度的白噪声组成，从而消除了模态混叠的困境[24-25]。
EEMD将信号x(t)（ICG信号、ECG信号或加速度信号）与白噪声n(t)相加，生成新的信号y(t)，再进一步分解为IMF，具体如式(2)和式(3)所示：
|
|
其中，zi(t)为x(t)信号分解的第i个IMF，K为得到的IMF总数，α为白噪声的幅值，rn(t)为EEMD分解后的残余信号。
由于ICG信号主频率范围为［0.7, 6.0］ Hz，本研究只选择主频率在这个范围内的IMF。为了保证后续CCA处理的性能，要求选择所选IMF的通道数至少为5个，如果实际符合要求的IMF少于5个，保留主频率邻近的IMF作为补充。
将所有被选中的IMF做快速傅里叶变换，并取其幅度部分分别记为：s1, ICG, ···, sN, ICG、 s1, ECG, ···, sN, ECG、s1, ACC, ···, sN, ACC，其中N为被选中的IMF个数，N ≤ K。定义ICG、ECG以及加速度信号的特征矩阵，SICG = ［s1, ICG, ···, sN, ICG］，SECG = ［s1, ECG, ···, sN, ECG］，SACC =［s1, ACC, ···, sN, ACC］。SICG、SECG以及SACC是频域CCA处理的输入信号，两步频域CCA主要是为了计算SICG、SECG以及SACC之间相关性，再根据相关性系数进一步处理带有运动伪迹的ICG信号。
假定CCA的典型相关变量和可以用目标信号Sa和Sb来表示：、 ，其中，w1和w2为权值向量，μa和νb为典型相关变量[23]，最终CCA的目标函数定义如式(4)所示：
|
其中，r为典型相关系数向量，和分别为Sa和Sb的自协方差矩阵，为Sa和Sb的交叉协方差矩阵。假设得到N个典型相关系数，用向量r来代表典型相关系数向量，r = ［r1, r2, ···, rk, ···, rN］；其中，0 ≤ ri ≤ 1。A和B的N个典型相关变量分别记录在矩阵Ua和Vb，相应的权值向量分别存储在矩阵Ca和Db中。
对于第一步频域CCA，计算ICG和加速度信号分解的IMF之间的频域相关性系数，以及ECG和加速度信号分解的IMF之间的频域相关性系数。这里认为相关性系数较大的部分是运动引起的噪声，因此将相关性系数高（ri > 0.8）的成分对应的典型相关变量置0，来抑制运动伪迹。
对于第二步频域CCA，计算ICG和ECG信号分解的IMF之间的频域相关性系数。此时，认为相关性系数大的部分是两个信号的主频成分，因此将相关性低（ri < 0.9）的成分对应的典型相关变量置0，来进一步去除运动伪迹。
2.1.3 重构ICG信号
重构ICG信号， ICG频谱可以通过保留的权值系数矩阵CICG和对应的部分分量被置为0的典型相关矩阵UICG进行重构，如式(5)所示：
|
其中，HICG = ［h1, ICG, ···, hN, ICG］是ICG信号重构的频谱，hi, ICG是ICG信号第i个IMF去噪后的频谱向量。
因此，ICG信号的第i个IMF去噪后的时域信号zi, ICG就可以利用频谱向量hi, ICG信息和对应的相位信息，再通过傅里叶逆变换计算得出，如式(6)所示：
|
其中，zi, ICG为第i个IMF去噪后的时域信号，M为傅里叶逆变换时的点数，π为圆周率，k为频域信号的频率索引，j为虚数单位，n为时域信号的时间索引。最后，将所有得到的时域IMF相加获得最终重构的ICG信号，如式（7）所示：
|
其中，yICG为最终重构的ICG信号，N为IMF的个数，n为时域信号的时间索引。
2.2 信号质量的检验指标
为了评估ICG信号去噪的有效性，需要对信号质量做进一步的评估，这是因为ICG的信号质量影响特征点检测的准确性，进而影响相应生理参数的计算准确性。由于缺乏评估ICG信号的量化指标，本文给出了通过ICG信号得到的相关生理参数的合理范围，如表1所示，并通过该范围限制来判定信号在去噪后的效果。一个有效的心拍应满足以下条件：① 射血前期（pre-ejection period，PEP）在50～160 ms的范围内；② 左室射血时间（left ventricular ejection time，LVET）在160～380 ms的范围内；③ 阻抗微分最大值在0.4～3.0 Ω/s的范围内；④ HR在30～180 beat/min的范围内[26]。
表1
ICG信号特征参数的正常值范围
Table1.
Normal range of ICG signal characteristic parameters
为了进一步验证信号质量，本文增加了更多细化的评价指标，主要包括：① 有效节拍（valid beat，VB）；② 节拍贡献因子（beat contribution factor，BCF），其计算方式为VB数目与总节拍数的比值，BCF越大说明信号质量越好；③ CC，为ICG信号中两个相邻C点的时间间隔，反映的是心动周期；CC计算出的心动周期与实际心动周期差值越小，信号质量越高；④ PEP变异性和LVET变异性，即在连续ICG数据的片段内（10 s），确认每个片段内测量的PEP和LVET的差异是否较小，亦即方差是否较小。
这些评价指标的测量依赖对ICG信号C点、B点以及X点等特征点的识别。ICG的C点是最大射血速度的时刻，是ICG在一个心动周期的幅度值最大的点[27-28]。ICG的B点是主动脉瓣的打开时刻，是幅值为C点幅值的0.15倍的特征点[29-31]。ICG的X点是主动脉瓣的关闭时刻，是ICG在一个心动周期的幅度值最小的点[27]。ICG信号的特征点检测方法可以总结如下：
（1）搜索ECG信号，识别特征点Q、R、T。
（2）扫描对应时刻的ICG信号，从ECG的R点对应时刻开始的三分之一个心动周期内，搜索ICG信号幅度的最大值，这个最大值对应特征点记为C点。
（3）扫描对应时刻的ICG信号，从ECG的R点对应时刻开始到ICG的C点对应时刻，找到幅值为C点幅值的0.15倍的特征点，并将其记为B点。
（4）扫描对应时刻的ICG信号，从ECG的T点对应时刻开始的三分之一心动周期内，将幅度最小值的特征点标记为X点。
（5）重复上述步骤，直到ICG信号结束。
3 结果与分析
3.1 结果
如图4所示，本文展示了ICG信号去噪效果和特征点识别检测，同时对ICG信号进行归一化，以便清晰地展示原始ICG和去噪后的ICG之间的差异。ICG信号的不同特征点用不同的颜色和标记来表示，可以清晰发现去噪后的信号获得了更精确的特征点。如表2所示，为30名受试者的测量结果。其中，静息测量场景共1 068次心拍；去噪后，ICG的VB数量从976次增加到1 055次，ICG的BCF也从原来的91.4%上升到98.7%。在运动测量场景中也有类似的改善，特别是ICG的BCF由78.4%提高到97.5%。总体而言，去噪前ICG信号的平均BCF约为80.1%，去噪后ICG的平均BCF约为97.7%。
图4
ICG信号去噪效果以及特征点识别
Figure4.
Denoising effect of ICG signal and feature point detection
表2
静止和运动两种状态下去噪效果的结果
Table2.
Results of the denoising effect in rest and motion state
由于ECG信号的R波特征明显，即使在运动伪迹干扰下也较为容易识别和测量，可以作为ICG信号的一个参考。如图5所示，去噪前估算出的心动周期与实际心动周期之间存在较大的误差值，尤其是在运动状态下。相比之下，去噪后估算出的心动周期和实际心动周期之间误差值明显减小，同时估计值与真实值的差值落在95%一致性区间内的数量更多。综上所述，在整个研究场景中（包括静息状态和运动状态），与去噪前相比，去噪后估算出的心动周期与实际心动周期更加一致，间接说明了ICG信号质量在去噪后得到提高。
图5
不同状态下心动周期的真实值与估算值的一致性
Figure5.
The consistence between the true value and the estimated value of the cardiac cycle in different states
如图6所示，对去噪后的ICG信号切片分组后，每个切片内数据获取的多个LVET值的方差均值明显低于去噪前。对于PEP值也有类似的结果，去噪后每个信号切片测量的PEP值的方差明显降低，间接说明去噪后ICG信号质量的提高。
图6
ICG信号在去噪前后，PEP和LVET参数的标准差的变化
Figure6.
Changes in the standard deviation of PEP and LVET parameters detected before and after ICG signal denoising
3.2 消融实验
上述实验结果验证了两步频域EEMD-CCA去噪方法对带有运动伪迹的ICG信号去噪的有效性。由于提出的去噪方法需要两步频域EEMD-CCA，本文通过消融实验来验证方法中每一步的必要性。第一步频域EEMD-CCA消除了噪声引起的主要频率成分，而第二步频域EEMD-CCA获取生理信号频谱中最相关的部分，进一步去除了运动产生的伪迹。如图7所示，单独使用第一步频域EEMD-CCA的去噪效果低于使用两步频域EEMD-CCA。如图8所示，单独使用第二步频域EEMD-CCA去噪效果也低于使用两步频域EEMD-CCA。因此，两步频域EEMD-CCA比仅仅使用其中一步频域EEMD-CCA在提高ICG质量方面更有效，验证了去噪方法的有效性。
图7
单独使用第一步频域EEMD-CCA与使用两步频域EEMD-CCA方法的去噪结果对比
Figure7.
Comparison of denoising results by using only the first step EEMD-CCA and the two-step EEMD-CCA method
图8
单独使用第二步频域EEMD-CCA方法与使用两步频域EEMD-CCA方法的去噪结果对比
Figure8.
Comparison of denoising results by using only the second step EEMD-CCA and the two-step EEMD-CCA method
4 总结
由于ICG信号的有效去噪有利于进一步的生理参数检测和心血管疾病诊断，故本文介绍了一种两步频域EEMD-CCA的方法对ICG信号进行去噪。本文方法在静息和运动两种状态下对30名受试者进行了评估。研究结果表明，去除运动伪迹后，ICG信号质量明显提高。
本文研究验证了所提去噪方法在静息和运动状态场景下的有效性，后续还需要进一步检验去噪方法在其他运动场景中的表现，如跑步、打球等运动幅度更大的场景下的有效性；同时，对ICG信号去噪只是心功能监测的一部分，后续还需要结合胸阻抗模型进一步验证相关心功能参数的准确性，如是否有利于对SV、CO计算准确度的提高。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：解尧负责软硬件平台搭建以及论文撰写、修订；杨东、余洪龙负责数据记录与分析；解启莲负责研究计划安排以及论文审阅修订。
伦理声明：本研究通过了安徽医科大学附属第四医院伦理委员会的审批（批号：PJ2021-006-01）。
