心冲击信号（ballistocardiogram，BCG）是表示心脏机械活动的间接信号，通过侦测血液泵入动脉引起的振动来获得. 随着传感器技术的发展，BCG信号^[1]可以在传感器不直接接触被测者的状态下进行采集. 该特性相较于心电图（electrocardiogram，ECG）信号采集有很大的区别，ECG通常需要将电极片贴在人的胸口、手臂、腿部等位置，被测者的异物感较强. 在BCG信号采集的生理床垫中，传感器位于心脏和肺部的正下方，被测者无异物感，舒适度高. 心脏跳动时产生的冲击力传导至柔性传感器，柔性传感器的电阻随着冲击力产生对应的变化，采集电路将电阻的变化转变成电压的变化. BCG信号的采集本质上是对压力信号的处理，可以认为是非平稳性信号^[2]. 从生理床垫传入的信号叠加了呼吸信号、体动伪影、外界干扰等信号. 其中，呼吸信号的电信号波幅值一般是心冲击信号的3~8倍^[3]，所有的生理信息通过生理床垫都转变成一维的压力信号^[4]，如何从这个一维模拟信号中去除噪声分量并提取有用生理信息成为研究的关键.

经验模态分解（empirical mode decomposition，EMD）方法由Huang等^[5]提出，作为1种传统的自适应分解信号的方法，被广泛应用于生物信号处理领域. 分解得到的分量称为固有模态函数（intrinsic mod functions，IMFs）.EMD的缺点是各个模态之间存在不同程度的模态混叠与端点效应，特别是在分离呼吸信号与心冲击信号这2种频率相近的信号时，问题尤为突出. 为了解决EMD模态混叠问题，Wu等^[6-7]改进了EMD算法，提出集合经验模态分解（ensemble empirical mode decomposition，EEMD）算法. 该方法取得了成功并大量应用于医疗^[8-9]、地质勘测^[10]、机械控制^[11]领域. EEMD通过加入白噪声，改变了信号极值点的分布，在一定程度上抑制了模态混叠现象. 在分解的过程中会有大量的伪分量混入，导致重构误差加大. Yeh等^[12]提出互补集合经验模态分解（complementary ensemble empirical mode decomposition，CEEMD）算法，通过多次加入符号相反的白噪声改善了该问题，但由于集成次数过多，导致计算效率低下、应用场景受到限制. 受到EMD的启发，Smith^[13]提出局部均值分解（local mean decomposition，LMD）的分解方法，但与EMD、EEMD、CEEMD方法一样，存在模态混叠、边界效应、欠包络等问题. 为了克服以上问题，Dragomiretskiy等^[14]提出变分模态分解（variational mode decomposition，VMD）的自适应分解方法. 该方法相对于EMD方法具有稳固的数学理论基础，VMD方法相比CEEMD、LMD等经典自适应方法具有良好的鲁棒性、收敛速度快的优点. 经分析大量的仿真信号发现，VMD算法对于低频信号的分离效果优于高频信号，此外VMD的分解精度受惩罚因子α和分解层数k的影响^[15]. 确定合适的k，才能避免模态混叠现象. 在现阶段，VMD方法用于分解生理信号的研究较少，属于研究的初始阶段.

在生理信号重构法中，目前常用的有两大类. 1）去噪类：利用小波变换选择基函数去除噪声^[16]或者是利用时序平移构造多维输入进行独立变量分析^[8]；2）提取类：在所有的IMF分量中，直接提取有效分量进行合成^[9]. 小波降噪方法虽然好，但会损失部分原始信号的信号特征且需要选择基函数；独立变量分析的联合去噪方法需要根据IMF选择信号与噪声的分界线来引入虚拟噪声通道，由于噪声与有效信号的分界线不固定. 分界线的判定会直接影响重构结果；直接提取法多采用经验所得的固定阈值^[4]，当环境变化时无法适用.

本文基于变分模态分解算法，分解数字去趋势（detrended fluctuation analysis，DFA）后的生理信号，改进能量提取法. 通过自适应筛选有效IMF相关分量，重构心冲击信号与呼吸信号，与EMD和CEEMD法进行重构效果对比.

BCG信号的典型波形如图1所示. 图中，H、I、J、K、L、M、N、O为单个心冲击信号周期内的8个极值点；其中H到K点为最主要的特征，反映了心脏收缩后血液击中主动脉和大血管的压力；L、M、N、O表示心脏舒张的过程.

生理信号分解重构流程如图2所示，按照先后顺序包括生理信号采集、硬件低通降噪、软件去趋势、生理信号VMD分解、自适应筛选有效分量、呼吸信号和心冲击信号的重构，共6个部分. 信号经ADC采样后由上位机Matlab中进行算法处理.

去趋势采用最优二乘拟合求得数据趋势线，用原始数据减去趋势线，达到去除原始数据基线漂移的目的. 处理后的数据位于X轴附近，方便了后期的数据处理. 具体步骤如下.

1）将N的时间序列x（t），t=1，2，···，N. 计算累计离差值并转换为序列，其中 $\overline x $为时间序列的平均值.

(1) $ S(t) = \sum\limits_{i = 1}^t {[x(i) - \bar x]}, $

(2) $ \overline x = \frac{1}{N}\sum\limits_{t = 1}^N {} x(t). $

2）将S_n（t）按照间距L划分为不重叠的n个相同长度区间，其中L为每一段区间长度，即时间尺度n为窗口数总和.

3）对每一段区间进行最小二乘线性拟合，得到每一段区间的变化趋势S_n（t）.

4）对S（t）减去每个独立区间的变化趋势，得到去趋势后的D（t）.

VMD通过迭代搜寻模型的最优解来确定每个分量的频率的中心和带宽，因此可以自适应地实现信号分解，将各分量分离成具有稀疏特性的独立分量. 在VMD算法中，拟定每个信号分量都有中心频率和带宽，所以变分问题即为寻求k个模态且使得模态带宽之和最小. VMD问题主要是求解下式的问题：

(3) $\left. \begin{aligned} &\mathop {\min }\limits_{{{u_k}},{{\omega _k}}}\; {\sum\limits_{k=1}^K {\left| {{\partial _t}\left[ {\left( {\delta (t) + \frac{{\rm{j}}}{{\text{π} t}}} \right)*{u_k}(t)} \right]{{\rm{exp}}\;{(- {\rm{j}}{\omega _k}t)}}} \right|} }; \\ &{\rm{s}}.{\rm{t}}{\rm{.}}\;\sum\limits_{k=1}^K {{u_k}}(t)= f(t). \end{aligned} \right\}$

式中： ${u_k}$为信号的第k个分解模态， ${\omega _k}$为每个模态信号的第k个中心频率.

为了求解上述变分问题，引入二次惩罚因子 $\alpha $与拉格朗日惩罚算子 $\lambda \left( t \right)$构造拉格朗日增广函数：

(4) $ \begin{split} L\left( {{u_k}{{,}}{\omega_k}{,}\lambda } \right): = & \alpha \sum\limits_{k=1}^K {\left| {{\partial _{{t}}}\left[ {\left( {\delta (t) + \frac{{\rm{j}}}{{\text{π} t}}} \right)*{u_k}(t)} \right]{{\rm{exp}}\;{ (- {\rm{j}}{\omega _k}t)}}} \right|} + \\ &\left| {f(t) - \sum\limits_{k=1}^K {{u_k}(t)} } \right| + \left\langle {\lambda (t),f(t) - \sum\limits_{k=1}^K {{u_k}(t)} } \right\rangle . \end{split} $

二次惩罚因子可以在独立同分布高斯噪声的影响下，保证了信号重构后的准确度. 拉格朗日乘数保证了过程的严格性.

采用乘法算子交替方向法（alternate direction method of multipliers，ADMM）交替更新从n步到n+1步迭代的 ${\widehat u_k^{n + 1}}$、 ${\omega _k^{n + 1}}$、 ${\widehat\lambda ^{n + 1}}$，将 ${\widehat u_k^{n + 1}}$、 ${\omega _k^{n + 1}}$的取值问题转化到频域，求得拉格朗日表达式的鞍点. 收敛条件为

(5) $\sum\limits_{k=1}^K {\frac{{\left| {\widehat u_k^{n + 1}(\omega)- \widehat u_k^n} (\omega)\right|}}{{\left| {\widehat u_k^n}(\omega) \right|}}} < \varepsilon .$

利用傅里叶变换，将迭代后的结果转换到每个模态的频域上进行更新，可得

(6) $\widehat u_k^{n + 1}(\omega ) = \frac{{\widehat f(\omega ) - \displaystyle\sum\limits_{i \ne k} {\widehat u_i^{n}(\omega ) + {{{\widehat \lambda}^{n} (\omega )}}/{2}} }}{{1 + 2\alpha {{(\omega - {\omega _k})}^2}}},$

(7) $\omega _k^{n + 1} = \frac{{\int_0^\infty \omega {{\left| {\widehat u_k^n(\omega )} \right|}^2}{\rm{d}}\omega }}{{\int_0^\infty {} {{\left| {\widehat u_k^n(\omega )} \right|}^2}{\rm{d}}\omega }},$

(8) $ \widehat \lambda _{}^{n + 1}(\omega ) = \widehat \lambda _{}^n(\omega ) + \tau \left[ {\widehat y_{}^n(\omega ) - \sum\limits_n {\widehat u_k^{n + 1}(\omega )} } \right]. $

式中： $\widehat f(\omega ) $ 为 $f(t) $ 的频域表示， $\widehat u_k^{n + 1}$为对当前剩余信号 $\widehat f(\omega ) - \displaystyle\sum\nolimits_{i \ne k} {\widehat u_i^{n}(\omega ) + {{{\widehat \lambda}^{n} (\omega )}}/{2}}$的维纳滤波； $\omega _k^{n + 1}$为当前模态功率谱的中心； $\widehat \lambda _{}^{n + 1}(\omega )$为对应的傅里叶变换. VMD算法的流程如下.

1）初始化 $\widehat u_k^1{\text{、}} \omega _k^1{\text{、}}\widehat \lambda ^1$，通过 Hilbert 变换得到每个模态函数的解析信号，获得信号的单边频率谱.

2）在第n+1步时，每一个模态函数都围绕各自估计的中心频率，通过加入指数项调整各自估计的中心频率，将每一个模态的频谱移动到基带上. 利用式（6）~（8），从1到K循环更新 ${\widehat u_k^{n+1}}$、 ${\omega_k^{n+1}}$、 $\widehat\lambda^{n+1}$.

3）当满足收敛条件（5）的约束时，输出K个窄带模态分量.

生理信号经过柔性压力传感器转换为数字信号时，包含了呼吸信号、心冲击信号、体动伪影信号、外界干扰等各种不同频段信号叠加的一维模拟信号. 在静息时，正常的呼吸频率主要集中在0.1~0.5 Hz，心冲击频率集中在1.0~3.0 Hz. 每个单位区间的生理信号都可被VMD分解为 ${\rm{IM}}{{\rm{F}}_1}\left( t \right), {\rm{IM}}{{\rm{F}}_2}\left( t \right),\cdots \;{\rm{IM}}{{\rm{F}}_K}\left( t \right)$，总共K个模态分量，呼吸与心冲击信号都至少由一个IMF分量构成. S_b（t）和S_h（t）分别表示重构后的呼吸信号和心冲击信号：

(9) ${S_{\rm{b}}}(t) = \sum\limits_{n = {n_{\rm{a}}}}^{{n_{_{\rm{b}}}}} {{\rm{IMF}}(t)} ,$

(10) ${S_{\rm{h}}}(t) = \sum\limits_{n = {n_{\rm{c}}}}^{{n_{\rm{d}}}} {{\rm{IMF}}(t)} .$

式中： ${n_{\rm{a}}} \leqslant {n_{\rm{b}}} \leqslant {n_{\rm{c}}} \leqslant {n_{\rm{d}}}$，等号不同时成立.

为了能够准确重构出呼吸信号和心冲击信号，需要确定n_a、n_b、n_c、n_d. 通过快速傅里叶变换分别计算每个IMF分量的能量 $E\left( i \right)$、每个IMF分量在呼吸信号频带0.1~0.5 Hz上的能量 ${E_{\rm{b}}}\left( i \right)$和心冲击信号频带1.0~3.0 Hz上的能量 ${E_{\rm{h}}}\left( i \right)$，判断是否满足以下条件：

(11) ${{{E_{\rm{b}}}(i)}}/{{E(i)}} \geqslant {\theta _{\rm{b}}},$

(12) ${{{E_{\rm{h}}}(i)}}/{{E(i)}} \geqslant {\theta _{\rm{h}}}.$

式中： ${\theta _{\rm{b}}}$和 ${\theta _{\rm{h}}}$分别为呼吸信号和心冲击信号的阈值. 为了提高重构生理信号的鲁棒性，采用自适应阈值的方法确定 ${\theta _{\rm{b}}}$和 ${\theta _{\rm{h}}}$. 具体流程如下.

1）分别找到呼吸信号和心冲击信号频带中能量占比最大的 $\max {E_{\rm{b}}}\left( i \right)$和 $\max{E_{\rm{h}}} \left( i \right)$.

2）筛选出 ${m_1}$个 ${E_{\rm{b}}}\left( i \right) > 0.1{E_{\rm{b}}}\max \left( i \right)$的呼吸信号有效阈值参考量 ${E_{{\rm{b\_ref}}}}(i)$、 ${m_2}$个 ${E_{\rm{h}}}\left( i \right) > 0.1{E_{\rm{h}}}\max \left( i \right)$的心冲击信号有效阈值参考量 ${E_{{\rm{h\_ref}}}}(i)$.

3）对求得的有效阈值参考量取算数平均值，求得当前的呼吸信号阈值 ${\theta _{\rm{b}}} = \displaystyle\sum\nolimits_{t = 1}^{{m_1}} {} {E_{{\rm{b\_ref}}}}(i)/{m_1}$，心冲击信号阈值 ${\theta _{\rm{h}}} =\displaystyle\sum\nolimits_{i = 1}^{{m_2}} {} {E_{{\rm{h\_ref}}}}(i)/{m_2}$.

取得阈值后，如果分解得到分量满足式（11），则认为该分量是组成呼吸信号的有效成分. 若分解得到的分量满足式（12），则认为该分量是组成呼吸信号的有效成分. 累加筛选出的IMF分量，得到呼吸信号和心冲击信号.

采用的24位ADC芯片为Texas Instruments公司的ads1220，采样频率设定为90 Hz，共采集4 096个点，共45.51 s. 柔性电阻型压力传感器采用Interlink Electronics公司的FSR408，压力为1~100 N. 主控采用STM32F407. 滤波采用巴特沃斯反馈型2阶低通滤波，截止频率为100 Hz，实物如图3所示. 被测者年龄25岁，男性，身高为180 cm，体重为75 kg. 传感器位于心脏的正下方，通过在传感器上以每5 cm间隔附着填充物，增加传感器的信号捕捉面积. 考虑到实际常出现的应用场景，采用一层薄床单与一件薄衣物厚度和约为4 mm. 睡姿采用最常见的仰卧姿势.

由于柔性压力传感器具有一定的延展性，在实际测量的过程中，当被测者躺在上面时必定会产生信号基线的漂移，如图4所示. 图中，V为信号的振幅. 信号需要经过3~10 s的时间，才会达到相对稳定的状态. 在采集数据上体现为一定程度的趋势波动. 若直接进行数据分析会引入误差，则需要对采集到的生理信号进行预处理，达到去除原始数据在稳定过程中的偏移分量、减小误差的目的.

DFA的过程本质上是将信号中的低频分量及直流量分离的过程. 为了防止DFA对正常数据结果造成影响，L的取值应大于最大呼吸间隔10 s对应的采样点900. L过大，则影响DFA的效果. 取L为1 024.

为了对比VMD算法在生理信号分解上的效果，使用EMD算法和CEEMD算法分解同一段生理信号并进行对比. 由分解结果可以看出，原始信号经过了EMD分解，得到9个从高频到低频的固有模态分量和1个趋势项，总共10个分量. 由于EMD的特性，分解时不需要预先设定参数，得到的IMF分量个数完全由信号本身特性决定. 为了对比效果，设定CEEMD与VMD分解层数为10. CEEMD采用默认设置，信噪比为0.2，迭代20次.

VMD算法中，惩罚因子α本质上为模态的中心约束强度，若分离的目标信号频率越接近，则需要的约束强度越大，即需要的通带带宽越小.α与通带带宽成反比. 经过多次试验可知，从0~40 000以1 000为步进值，在α小于11 000与大于29 000时较易产生混叠现象. 在层数k=10的前提下设定α为20 000. 如图5（a）~（c）所示分别为利用3种方法分解后的3组IMF分量.

EMD与CEEMD 的分解结果为IMF₁~IMF₁₀，分解顺序从高频开始到低频结束. VMD分解与前两者不同，分解顺序从低频开始到高频结束. 通过对比EMD与CEEMD的分解结果可知，EMD在第4个IMF与第5个IMF、第6个IMF与第7个IMF之间存在较明显的频谱混叠. CEEMD分解中得益于辅助噪声的叠加，混叠现象有所缓解，但对呼吸信号和心冲击信号的重构产生了影响. VMD在分解效果上优于前两者，从图5（c）的分解结果来看，VMD分解所得的IMF中，大多都集中在目标生理信号的0.1~3.0 Hz内，且没有大面积的频谱重叠.

如图6所示为利用3种分解方法分解出的分量通过HHT变换得到的Hilbert-Huang谱^[17-19]. 在理想情况下，每个IMF分量之间都应完全正交. 从图6可以看出，利用3种方法均能将信号分解成若干频率成分. EMD和CEEMD分解中极值点分布不均匀，分解过程中的包络为异常包络和正常包络的混合，产生了不同程度的模态混叠和端点效应. 在HHT后的频谱表现为相邻谱线相交和同一条谱线的分布较分散. 在生理信号分离实验中，EMD分解和CEEMD分解后在频率大于0.5 Hz的IMF分量上出现了明显的谱线分散. 在VMD分解中，各个IMF分量都相对独立，未出现大量重叠. VMD分解在一定程度上克服EMD和CEEMD分解过程中存在的模态频谱混叠和边界效应的问题.

阈值 ${\theta _{\rm{b}}}$和 ${\theta _{\rm{h}}}$决定了重构的生理信号的准确度.图7和表1给出利用3种方法分解所得IMF各分量在呼吸信号频带和心冲击信号频带内能量的百分比P_b与P_h. 根据自适应确定阈值的方法，确定重构分量.

EMD法分解中，呼吸信号P_bmax为IMF₆的0.820 7，心冲击信号P_hmax为IMF₄的0.577 7，在滤除远小于最大能量的无效阈值后，呼吸信号中有效的阈值参考量为IMF₅、IMF₆、IMF₇、IMF₈. 心冲击信号中有效的阈值参考量为IMF₂、IMF₃、IMF₄、IMF₅、IMF₇、IMF₈. 分别对上述2组IMF能量求均值，求得本次EMD分解后用于重构信号的阈值 ${\theta _{\rm{b}}}$=0.414 6， ${\theta _{\rm{h}}}$=0.277 2. 当IMF分量中存在的生理信号能量占比达到阈值时，把该IMF分量确定为待重构分量. 分别将阈值代入式（11）、（12），可得EMD法中用于重构呼吸信号的分量有IMF₆与IMF₇ 2个分量，用于重构心冲击信号分量的有IMF₃、IMF₄、IMF₅ 3个分量.

CEEMD法分解中，呼吸信号P_bmax为IMF₅的0.901 7，心冲击信号P_hmax为IMF₃的0.442 4. 呼吸信号中有效的阈值参考量为IMF₅、IMF₆、IMF₇. 心冲击信号中有效的阈值参考量为IMF₂、IMF₃、IMF₄、IMF₅、IMF₆、IMF₇. 求得本次CEEMD分解后用于重构信号的阈值 ${\theta _{\rm{b}}}$=0.547 4， ${\theta _{\rm{h}}}$=0.252 1. 可得用于重构呼吸信号的分量有IMF₅与IMF₆ 2个分量，用于重构心冲击信号分量的有IMF₃、IMF₄ 2个分量.

VMD法分解中，呼吸信号P_bmax为IMF₁的0.858 1，心冲击信号P_hmax为IMF₃的0.394 0，呼吸信号中的有效阈值参考量为IMF₁、IMF₂. 心冲击信号中有效阈值参考量为IMF₁、IMF₂、IMF₃、IMF₄、IMF₅、IMF₆. 求得本次VMD分解后用于重构信号的阈值 ${\theta _{\rm{b}}}$=0.554 9， ${\theta _{\rm{h}}}$=0.217 0. 可得用于重构呼吸信号的分量有IMF₁ 1个分量，用于重构心冲击信号分量的有IMF₃、IMF₄、IMF₅ 3个分量.

如图8（a）所示为利用3种分解方法重构后的心冲击信号与对应ECG信号效果图. ECG信号的采集使用PC-80B快速心电检测仪，采用V5胸前导联接线法，通过时间戳同步ECG信号与BCG信号. 基于VMD重构后的心冲击信号可以看到H、I、J、K共4个极值点. 在基于EMD及CEEMD的分解方法中，由于频谱的混叠，心冲击信号的重构混入了3 Hz以上频率的干扰，在信号上表现出了不规则的波动，只有J极值点较明显. 3种分解方式在实际的测试中，对比心冲击信号典型波形发现，信号的L、M、N、O极值点不能明显地还原. 因为在心脏的舒张过程中柔性压力传感器与人体不是完美的贴合，在传感器与皮肤之间的衣物与床单会将部分微弱的生理信号吸收.

如图8（b）所示为3种方法的呼吸信号重构波形与原始信号对比. 呼吸频率的标准值采用人工标定计算，当开始吸气时为1个呼吸周期的开始. 基于EMD自适应分解方法重构出的呼吸信号波形在30 s处出现了信号缺失. 基于CEEMD法中的呼吸信号因为低频分量的混入使重构波形出现了上下起伏. 基于VMD法重构的呼吸信号无信号缺失与上下波动，与原始信号中的呼吸波形最相似.

为了检验本文重构方法的准确性，采集80组原始数据，参数与上文设置一致，依次用EMD、CEEMD、VMD 3种方法分解重构，得到心冲击信号和呼吸信号. 经FFT计算得到的心率和呼吸率作为测量值数据. 使用快速心电检测仪与人工计数的方法，获得心率和呼吸率的标准值. 使用Bland-Altman分析法^[20-21]判断分解重构的一致性，分析结果如图9所示. 图中，X为标准值与当前测量值的平均值，Y为标准值与当前测量值的差值，每个点表示一次测量. 在X轴下方的直线表示80次实验的差值的平均值 $\bar d $，在X轴上、下的2条虚线表示95%一致性标准界限的上、下限 $\bar d $±1.96S_d（其中S_d为标准差）. 当表示差值均数的实线越接近X轴，一致性界限外的数据越少时，说明2种测量方法的一致程度越高.

利用EMD法得到的呼吸信号测量值与标准值差值的均数为−0.6，一致性界限为[−3.5，2.3]，共有90%的点在95%的置信区间内. 利用EMD法得到的心冲击信号测量值与标准值差值的均数为−0.3，一致性界限为[−4.5，3.9]，有91.25%的点在95%的置信区间内.

利用CEEMD法得到的呼吸信号测量值与标准值差值的均数为−0.2，一致性界限为[−3.2，2.7]，共有91.25%的点在95%的置信区间内. 利用CEEMD法得到的心冲击信号测量值与标准值差值的均数为−0.2，一致性界限为[−3.9，3.5]，有92.5%的点在95%的置信区间内.

利用VMD法得到的呼吸信号测量值与标准值差值的均数为−0.12，一致性界限为[−1.91，1.66]，共有92.5%的点在95%的置信区间内. 利用VMD法得到的心冲击信号测量值与标准值差值的均数为−0.2，一致性界限为[−3.7，3.4]，共有93.75%的点在95%的置信区间内. 以上数据结果表明，采用的通过VMD分解并重构呼吸信号与心冲击信号的方法相对于传统分解方法有较高的一致性.

（1）提出非接触式生理信号监测床垫.

（2）基于VMD的生理信号分解，有效改善了传统EMD与CEEMD在处理频段相近的生理信号时出现的频谱混叠问题.

（3）提出基于IMF能量分析自适应选取有效分量方法，提高了重构生理信号的准确性，比传统方法中固定阈值法的适应性更好.

（4）通过Bland-Altman分析法可知，基于VMD法比EMD与CEEMD法在呼吸与心冲击信号重构的应用中有更高的一致性.

本文将VMD应用于睡眠监测床垫中的生命体征信号提取，通过自适应方法筛选分解后的IMF分量在呼吸信号和心冲击信号中的有效分量，重构了BCG与呼吸信号，为非接触式睡眠监测的应用提供了新的思路. 虽然VMD算法在效果上优于EMD与CEEMD，但需要进一步的参数优化来缩短计算时间.