一种基于信息增益的心肺耦合关系分析方法与流程

1.本发明涉及心肺耦合关系分析的技术领域，尤其涉及一种基于信息增益的心肺耦合关系分析方法。


背景技术：

2.随着人们生活、工作的压力越来越大，引发的心血管疾病率也呈现出逐年持续上升的趋势，心脏方面的疾病变成了造成死亡的主要因素。很多国家都有数以万计的人患有不同程度的心血管方面的疾病，而在心脏病的诊断治疗中，使用最广泛最重要的检测技术是心电图，但单一的分析心电信号或者呼吸信号，分析结果的准确性和稳定性会受到序列非稳性的严重影响，如何实现心肺耦合关系分析，从多角度研究心血管疾病，成为当前研究的热门话题。
3.鉴于此，本发明提出一种基于信息增益的心肺耦合关系分析方法，通过对心电信号序列以及呼吸信号序列进行经验模态分解，对分解后的心电信号序列和呼吸信号序列进行信息增益计算，信息增益越小，则说明心肺耦合关系越强，从而实现心肺耦合关系分析。


技术实现要素：

4.本发明提供一种基于信息增益的心肺耦合关系分析方法，目的在于(1)实现心电信号序列以及呼吸信号序列的分解；(2)基于信息增益实现心肺耦合关系分析。
5.实现上述目的，本发明提供的一种基于信息增益的心肺耦合关系分析方法，包括以下步骤：
6.s1：采集人体的心电信号和呼吸信号，形成心电信号序列和呼吸信号序列，分别对采集得到的心电信号序列和呼吸信号序列进行预处理；
7.s2：利用变步长的自适用滤波算法对预处理后的心电信号序列和呼吸信号序列进行滤波处理；
8.s3：对滤波后的心电信号序列和呼吸信号序列进行经验模态分解；
9.s4：对分解后的心电信号序列和呼吸信号序列进行信息增益计算，信息增益越小，则说明心肺耦合关系越强，否则越低。
10.作为本发明的进一步改进方法：
11.所述s1步骤中采集人体的心电信号序列以及呼吸信号序列，包括：
12.将心电电极片固定在人体皮肤部位，心电电极片持续获取人体的心电信号，并将心电信号传输到心电传感器，心电传感器所接收的心电信号序列为d(t)；在心电传感器内部具有高通滤波器，高通滤波器允许心电传感器接收信号频率高于截止频率f
t
的心电信号，所述截止频率f
t
的计算公式为：
[0013][0014]
其中：
[0015]
r表示高通滤波器中电阻的电阻值；
[0016]
c表示高通滤波器中电容器的电容；
[0017]
将二极管作为温度传感器放置在人体鼻腔内，将由鼻腔内温度变化引起的二极管导通电压变化转换为交变电压信号，并对交变电压信号进行放大处理，得到放大后的呼吸信号序列u(t)，所述交变电压信号放大处理流程为：
[0018]
将采集到的交变电压信号u1(t)输入到信号运算放大器的同相输入端，则同相输入端的信号输出结果u2(t)为：
[0019]
u2(t)＝(r1 r2)u1(t)
[0020]
其中：
[0021]
r1,r2分别为信号运算放大器的同相电阻，r1 r2》1；
[0022]
将同相输入端的输出结果u2(t)输入到信号运算放大器的反相输入端，则反相输入端的信号输出结果为：
[0023]
u(t)＝u2(t)/r3[0024]
其中：
[0025]
r3为信号运算放大器的反相电阻，r3》1。
[0026]
所述s1步骤中对采集到的心电信号序列以及呼吸信号序列进行预处理，得到预处理后的心电信号序列和呼吸信号序列，包括：
[0027]
对采集到的信号x(t)进行膨胀腐蚀操作，得到降噪后的信号x
′
(t)，其中信号x(t)包括采集到的心电信号序列和呼吸信号序列，所述数学形态操作流程为：
[0028]
对采集到的信号x(t)进行g个点的膨胀腐蚀运算处理：
[0029][0030][0031][0032]
其中：
[0033]
表示信号的腐蚀运算，该运算使得x(t)信号序列在信号前段减少了g-1个信号点；
[0034]
表示信号的膨胀运算，该运算使得x(t)信号序列在信号后段减少了g-1个信号点；
[0035]
b(g)表示结构元素，g∈{0,1,
…
,g-1}，且g＜＜n；
[0036]
x(t)表示滤波处理后的信号，t∈{0,1,
…
,n-1}，n表示信号x(t)中的时刻总数；
[0037]
x
′
(t)表示降噪后的信号，包括降噪后的心电信号序列和呼吸信号序列。
[0038]
所述s2步骤中对预处理后的心电信号序列和呼吸信号序列进行自适用滤波处理，包括：
[0039]
利用变步长的自适用滤波算法对预处理后的心电信号序列和呼吸信号序列进行滤波处理，所述变步长的自适用滤波算法流程为：
[0040]
1)对于降噪处理后的信号序列x
′
(t)，令m∈t，m＝0,1,2,
…
，得到x
′
(m)，m的初始值为0，步长μ(0)＝0；
[0041]
2)输出滤波后的信号y(m)：
[0042]
y(m)＝x
′
(m)
t
ω(m)
[0043]
其中：
[0044]
ω(m)为滤波权值；
[0045]
t表示转置；
[0046]
3)计算信号序列误差：
[0047]
e(m)＝q(m)-ω(m)
t
x
′
(m)
[0048]
其中：
[0049]
q(m)为期望输出的滤波信号；
[0050]
4)更新步长：
[0051]
μ(e(m))＝αμ(e(m-1)) (1-α)a|e(m)|2exp(-b|e(m)|2)
[0052]
其中：
[0053]
α为平滑因子，α值越大步长变化越稳定，受误差e(m)的影响越小，将其设置为0.8；
[0054]
a,b为常数，将a设置为2，b设置为3；
[0055]
5)更新权值：
[0056]
ω(m 1)＝ω(m) μ(e(m))|e(m)|
v-2
e(m)x
′
(m)
[0057]
其中：
[0058]
v∈(0,2)，为平均范数；
[0059]
更新m＝m 1；
[0060]
6)重复步骤2)-5)，得到自适用滤波处理后的信号序列y(t)，其中y(t)为自适用滤波处理后的心电信号序列或呼吸信号序列。
[0061]
所述s3步骤中对自适用滤波后的心电信号序列和呼吸信号序列进行经验模态分解，包括：
[0062]
对自适用滤波后的心电信号序列和呼吸信号序列进行经验模态分解，所述经验模态分解算法流程为：
[0063]
1)将信号序列y(t)加入同长度的白噪声n(t)，得到信号y(t)：
[0064]
y(t)＝y(t) n(t)
[0065]
其中：
[0066]
y(t)为自适用滤波处理后的心电信号序列或呼吸信号序列；
[0067]
2)遍历出信号y(t)中所有的极大值点和极小值点，将所有的极大值点以及所有的极小值点分别用三次样条曲线连接，分别得到信号y(t)的上包络线以及下包络线，上下包络线的均值曲线即为m1(t)，则h1(t)＝y(t)-m1(t)；
[0068]
3)对h1(t)进行步骤2)操作，得到h2(t)；
[0069]
4)重复步骤3)，得到{h1(t),h2(t),
…
,h
10
(t)}，将所得到的h
10
(t)作为信号y(t)的第一组imf分量c1(t)，第一组imf分量的剩余分量r1(t)＝y(t)-c1(t)；
[0070]
5)判断剩余分量ri(t)中的i是否大于k，若i≤k，对ri(t)进行步骤2)-4)，即将ri(t)作为步骤2)中待遍历的信号，步骤4)中输出的剩余分量为r
i 1
(t)，并输出第i组imf分量，若i》k,则转到步骤6)；
[0071]
6)得到k-1组imf分量c2(t),c3(t),
…
,ck(t)以及剩余分量rk(t)＝r
k-1
(t)-ck(t)；
则将信号y(t)分解为下式：
[0072][0073]
将分解得到的imf分量{c1(t),c2(t),
…
,ci(t),
…
,ck(t)}进行快速傅里叶变换，分别得到呼吸信号和心电信号imf分量的频谱：
[0074]
ffi＝fft(ci(t))
[0075]
其中：
[0076]
ffi表示第i组imf分量的频谱；
[0077]ci
(t)表示信号y(t)的第i组imf分量；
[0078]
fft(
·
)表示快速傅里叶变换处理；
[0079]
对于心电信号序列，若imf分量频谱中0-60hz的频域能量达到整个频域能量的90％以上，则将该imf分量判断为心电信号的imf分量，否则将该imf分量判断为噪声信号的imf分量，并删除噪声信号的imf分量，得到心电信号序列的imf分量集合：
[0080]
q＝{c
1,1
(t),c
1,2
(t),
…
,c
1,q
(t)}
[0081]
其中：
[0082]
q表示心电信号序列中的imf分量数目；
[0083]
对于呼吸信号序列，若imf分量频谱中0-45hz的频域能量达到整个频域能量的90％以上，则将该imf分量判断为呼吸信号的imf分量，否则将该imf分量判断为噪声信号的imf分量，并删除噪声信号的imf分量，得到呼吸信号序列的imf分量集合：
[0084]
p＝{c
2,1
(t),c
2,2
(t),
…
,c
2,p
(t)}
[0085]
p表示呼吸信号序列中的imf分量数目。
[0086]
所述s4步骤中计算分解后心电信号序列和呼吸信号序列的信息增益，包括：
[0087]
计算心电信号序列分解后的imf分量与呼吸信号序列分解后的imf分量之间的信息增益h(p,q)：
[0088][0089]
其中：
[0090]
p
ij
表示呼吸信号在第i个imf分量c
2,i
(t)状态时，心电信号在第j个imf分量c
1,j
(t)状态的概率，在本发明一个具体实施例中，p
ij
表示c
2,i
(t)与c
1,j
(t)之间的余弦相似度。
[0091]
所述s4步骤中比较信息增益与预设阈值的大小，若信息增益小于预设阈值，则说明心肺具有较强的耦合关系，信息增益越小，心肺耦合关系越强，包括：
[0092]
预设耦合关系阈值对于人体心肺的信息增益值h(p,q)，若则说明心肺具有较强的耦合关系，其中信息增益越小，心肺耦合关系越强。
[0093]
相对于现有技术，本发明提出一种基于信息增益的心肺耦合关系分析方法，该技术具有以下优势：
[0094]
首先，本方案提出一种变步长的自适用滤波算法对预处理后的心电信号序列和呼吸信号序列进行滤波处理，所述变步长的自适用滤波算法流程为：对于降噪处理后的信号
序列x
′
(t)，令m∈t，m＝0,1,2,
…
，得到x
′
(m)，m的初始值为0，步长μ(0)＝0；输出滤波后的信号y(m)：
[0095]
y(m)＝x
′
(m)
t
ω(m)
[0096]
其中：ω(m)为滤波权值；
[0097]
t表示转置；
[0098]
计算信号序列误差：
[0099]
e(m)＝q(m)-ω(m)
t
x
′
(m)
[0100]
q(m)为期望输出的滤波信号；更新步长：
[0101]
μ(e(m))＝αμ(e(m-1)) (1-α)a|e(m)|2exp(-b|e(m)|2)
[0102]
其中：α为平滑因子，α值越大步长变化越稳定，受误差e(m)的影响越小，将其设置为0.8；a,b为常数，将a设置为2，b设置为3；更新权值：
[0103]
ω(m 1)＝ω(m) μ(e(m))|e(m)|
v-2
e(m)x
′
(m)
[0104]
v∈(0,2)，为平均范数；更新m＝m 1；重复上述步骤步骤，得到自适用滤波处理后的信号序列y(t)，其中y(t)为自适用滤波处理后的心电信号序列或呼吸信号序列。相较于传统方案，本方案所述步长μ(e(m))在滤波器受到脉冲噪声干扰时，若脉冲噪声较大，μ(e(m))值很小且趋近于0，μ(e(m))在当误差在附近时能取得较大值，以提高收敛速度；当误差e(m)接近原点时，μ(e(m))逐渐接近于0，以降低稳态误差，从而实现信号序列的自适用滤波处理。
[0105]
同时，本方案提出一种心肺耦合关系分析方法，通过对自适用滤波后的心电信号序列和呼吸信号序列进行经验模态分解，分解得到若干imf分量，将分解得到的imf分量{c1(t),c2(t),
…
,ci(t),
…
,ck(t)}进行快速傅里叶变换，分别得到呼吸信号和心电信号imf分量的频谱：
[0106]
ffi＝fft(ci(t))
[0107]
其中：ffi表示第i组imf分量的频谱；ci(t)表示信号y(t)的第i组imf分量；fft(
·
)表示快速傅里叶变换处理；对于心电信号序列，若imf分量频谱中0-60hz的频域能量达到整个频域能量的90％以上，则将该imf分量判断为心电信号的imf分量，否则将该imf分量判断为噪声信号的imf分量，并删除噪声信号的imf分量，得到心电信号序列的imf分量集合：
[0108]
q＝{c
1,1
(t),c
1,2
(t),
…
,c
1,q
(t)}
[0109]
其中：q表示心电信号序列中的imf分量数目；对于呼吸信号序列，若imf分量频谱中0-45hz的频域能量达到整个频域能量的90％以上，则将该imf分量判断为呼吸信号的imf分量，否则将该imf分量判断为噪声信号的imf分量，并删除噪声信号的imf分量，得到呼吸信号序列的imf分量集合：
[0110]
p＝{c
2,1
(t),c
2,2
(t),
…
,c
2,p
(t)}
[0111]
p表示心电信号序列中的imf分量数目，相较于传统方案，本方案分别对心电信号序列以及呼吸信号序列进行不同的经验模态分解，得到两种信号序列的imf分量，通过计算心电信号序列分解后的imf分量与呼吸信号序列分解后的imf分量之间的信息增益h(p,q)：
[0112]
[0113]
其中：p
ij
表示呼吸信号在第i个imf分量c
2,i
(t)状态时，心电信号在第j个imf分量c
1,j
(t)状态的概率。通过预设耦合关系阈值对于人体心肺的信息增益值h(p,q)，若则说明心肺具有较强的耦合关系，其中信息增益越小，心肺耦合关系越强，从而实现心肺耦合关系分析。
附图说明
[0114]
图1为本发明一实施例提供的一种基于信息增益的心肺耦合关系分析方法的流程示意图；
[0115]
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0116]
应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0117]
s1：采集人体的心电信号和呼吸信号，形成心电信号序列和呼吸信号序列，分别对采集得到的心电信号序列和呼吸信号序列进行预处理。
[0118]
所述s1步骤中采集人体的心电信号序列以及呼吸信号序列，包括：
[0119]
将心电电极片固定在人体皮肤部位，心电电极片持续获取人体的心电信号，并将心电信号传输到心电传感器，心电传感器所接收的心电信号序列为d(t)；在心电传感器内部具有高通滤波器，高通滤波器允许心电传感器接收信号频率高于截止频率f
t
的心电信号，所述截止频率f
t
的计算公式为：
[0120][0121]
其中：
[0122]
r表示高通滤波器中电阻的电阻值；
[0123]
c表示高通滤波器中电容器的电容；
[0124]
将二极管作为温度传感器放置在人体鼻腔内，将由鼻腔内温度变化引起的二极管导通电压变化转换为交变电压信号，并对交变电压信号进行放大处理，得到放大后的呼吸信号序列u(t)，所述交变电压信号放大处理流程为：
[0125]
将采集到的交变电压信号u1(t)输入到信号运算放大器的同相输入端，则同相输入端的信号输出结果u2(t)为：
[0126]
u2(t)＝(r1 r2)u1(t)
[0127]
其中：
[0128]
r1,r2分别为信号运算放大器的同相电阻，r1 r2》1；
[0129]
将同相输入端的输出结果u2(t)输入到信号运算放大器的反相输入端，则反相输入端的信号输出结果为：
[0130]
u(t)＝u2(t)/r3[0131]
其中：
[0132]
r3为信号运算放大器的反相电阻，r3》1。
[0133]
所述s1步骤中对采集到的心电信号序列以及呼吸信号序列进行预处理，得到预处理后的心电信号序列和呼吸信号序列，包括：
[0134]
对采集到的信号x(t)进行膨胀腐蚀操作，得到降噪后的信号x
′
(t)，其中信号x(t)包括采集到的心电信号序列和呼吸信号序列，所述数学形态操作流程为：
[0135]
对采集到的信号x(t)进行g个点的膨胀腐蚀运算处理：
[0136][0137][0138][0139]
其中：
[0140]
表示信号的腐蚀运算，该运算使得x(t)信号序列在信号前段减少了g-1个信号点；
[0141]
表示信号的膨胀运算，该运算使得x(t)信号序列在信号后段减少了g-1个信号点；
[0142]
b(g)表示结构元素，g∈{0,1,
…
,g-1}，且g＜＜n；
[0143]
x(t)表示滤波处理后的信号，t∈{0,1,
…
,n-1}，n表示信号x(t)中的时刻总数；
[0144]
x
′
(t)表示降噪后的信号，包括降噪后的心电信号序列和呼吸信号序列。
[0145]
s2：利用变步长的自适用滤波算法对预处理后的心电信号序列和呼吸信号序列进行滤波处理。
[0146]
所述s2步骤中对预处理后的心电信号序列和呼吸信号序列进行自适用滤波处理，包括：
[0147]
利用变步长的自适用滤波算法对预处理后的心电信号序列和呼吸信号序列进行滤波处理，所述变步长的自适用滤波算法流程为：
[0148]
1)对于降噪处理后的信号序列x
′
(t)，令m∈t，m＝0,1,2,
…
，得到x
′
(m)，m的初始值为0，步长μ(0)＝0；
[0149]
2)输出滤波后的信号y(m)：
[0150]
y(m)＝x
′
(m)
t
ω(m)
[0151]
其中：
[0152]
ω(m)为滤波权值；
[0153]
t表示转置；
[0154]
3)计算信号序列误差：
[0155]
e(m)＝q(m)-ω(m)
t
x
′
(m)
[0156]
其中：
[0157]
q(m)为期望输出的滤波信号；
[0158]
4)更新步长：
[0159]
μ(e(m))＝αμ(e(m-1)) (1-α)a|e(m)|2exp(-b|e(m)|2)
[0160]
其中：
[0161]
α为平滑因子，α值越大步长变化越稳定，受误差e(m)的影响越小，将其设置为0.8；
[0162]
a,b为常数，将a设置为2，b设置为3；
[0163]
5)更新权值：
[0164]
ω(m 1)＝ω(m) μ(e(m))|e(m)|
v-2
e(m)x
′
(m)
[0165]
其中：
[0166]
v∈(0,2)，为平均范数；
[0167]
更新m＝m 1；
[0168]
6)重复步骤2)-5)，得到自适用滤波处理后的信号序列y(t)，其中y(t)为自适用滤波处理后的心电信号序列或呼吸信号序列。
[0169]
s3：对滤波后的心电信号序列和呼吸信号序列进行经验模态分解。
[0170]
所述s3步骤中对自适用滤波后的心电信号序列和呼吸信号序列进行经验模态分解，包括：
[0171]
对自适用滤波后的心电信号序列和呼吸信号序列进行经验模态分解，所述经验模态分解算法流程为：
[0172]
1)将信号序列y(t)加入同长度的白噪声n(t)，得到信号y(t)：
[0173]
y(t)＝y(t) n(t)
[0174]
其中：
[0175]
y(t)为自适用滤波处理后的心电信号序列或呼吸信号序列；
[0176]
2)遍历出信号y(t)中所有的极大值点和极小值点，将所有的极大值点以及所有的极小值点分别用三次样条曲线连接，分别得到信号y(t)的上包络线以及下包络线，上下包络线的均值曲线即为m1(t)，则h1(t)＝y(t)-m1(t)；
[0177]
3)对h1(t)进行步骤2)操作，得到h2(t)；
[0178]
4)重复步骤3)，得到{h1(t),h2(t),
…
,h
10
(t)}，将所得到的h
10
(t)作为信号y(t)的第一组imf分量c1(t)，第一组imf分量的剩余分量r1(t)＝y(t)-c1(t)；
[0179]
5)判断剩余分量ri(t)中的i是否大于k，若i≤k，对ri(t)进行步骤2)-4)，即将ri(t)作为步骤2)中待遍历的信号，步骤4)中输出的剩余分量为r
i 1
(t)，并输出第i组imf分量，若i》k,则转到步骤6)；
[0180]
6)得到k-1组imf分量c2(t),c3(t),
…
,ck(t)以及剩余分量rk(t)＝r
k-1
(t)-ck(t)；则将信号y(t)分解为下式：
[0181][0182]
将分解得到的imf分量{c1(t),c2(t),
…
,ci(t),
…
,ck(t)}进行快速傅里叶变换，分别得到呼吸信号和心电信号imf分量的频谱：
[0183]
ffi＝fft(ci(t))
[0184]
其中：
[0185]
ffi表示第i组imf分量的频谱；
[0186]ci
(t)表示信号y(t)的第i组imf分量；
[0187]
fft(
·
)表示快速傅里叶变换处理；
[0188]
对于心电信号序列，若imf分量频谱中0-60hz的频域能量达到整个频域能量的90％以上，则将该imf分量判断为心电信号的imf分量，否则将该imf分量判断为噪声信号的imf分量，并删除噪声信号的imf分量，得到心电信号序列的imf分量集合：
[0189]
q＝{c
1,1
(t),c
1,2
(t),
…
,c
1,q
(t)}
[0190]
其中：
[0191]
q表示心电信号序列中的imf分量数目；
[0192]
对于呼吸信号序列，若imf分量频谱中0-45hz的频域能量达到整个频域能量的90％以上，则将该imf分量判断为呼吸信号的imf分量，否则将该imf分量判断为噪声信号的imf分量，并删除噪声信号的imf分量，得到呼吸信号序列的imf分量集合：
[0193]
p＝{c
2,1
(t),c
2,2
(t),
…
,c
2,p
(t)}
[0194]
p表示呼吸信号序列中的imf分量数目。
[0195]
s4：对分解后的心电信号序列和呼吸信号序列进行信息增益计算，信息增益越小，则说明心肺耦合关系越强，否则越低。
[0196]
所述s4步骤中计算分解后心电信号序列和呼吸信号序列的信息增益，包括：
[0197]
计算心电信号序列分解后的imf分量与呼吸信号序列分解后的imf分量之间的信息增益h(p,q)：
[0198][0199]
其中：
[0200]
p
ij
表示呼吸信号在第i个imf分量c
2,i
(t)状态时，心电信号在第j个imf分量c
1,j
(t)状态的概率，在本发明一个具体实施例中，p
ij
表示c
2,i
(t)与c
1,j
(t)之间的余弦相似度。
[0201]
所述s4步骤中比较信息增益与预设阈值的大小，若信息增益小于预设阈值，则说明心肺具有较强的耦合关系，信息增益越小，心肺耦合关系越强，包括：
[0202]
预设耦合关系阈值对于人体心肺的信息增益值h(p,q)，若则说明心肺具有较强的耦合关系，其中信息增益越小，心肺耦合关系越强。
[0203]
需要说明的是，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。并且本文中的术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、装置、物品或者方法不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、装置、物品或者方法所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括该要素的过程、装置、物品或者方法中还存在另外的相同要素。
[0204]
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在如上所述的一个存储介质(如rom/ram、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。
[0205]
以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。