基于二维空间的时间序列样本熵的计算方法及系统与流程

本发明涉及信号处理技术领域，特别涉及一种基于二维空间的时间序列样本熵的计算方法及系统。

背景技术：

在对时间序列信号进行处理时，时域的非线性分析方法是脑电信号处理的重要方面，包括相关维数、李雅普诺夫指数、小波熵、近似熵、样本熵等都是比较常用的方法。振幅和周期是表达波形振动的两个重要变量，现有的脑电可以依据波形的振动特性从时间序列中提取振幅和周期特征，并结合两种特征进行非线性分析。

样本熵是对振动时间序列的复杂度进行分析的一种指标，表征了时间序列中新模式的生成概率。以往的样本熵研究中，通常只将波形作为一维序列进行分析，而不能对波形的振幅-周期综合特性进行分析。

技术实现要素：

鉴于上述问题，提出了本发明以便提供一种克服上述问题或者至少部分地解决上述问题的一种基于二维空间的时间序列样本熵的计算方法及系统。

依据本发明的一个方面，提供了一种基于二维空间的时间序列样本熵的计算方法，所述方法包括：

S1：提取时间序列的极值点和极值点的时间位置，以生成极值序列P＝{p(i)}及时间位置序列T＝{t(i)}，其中，i＝1,2,3,…,n；

S2：生成序列Z₁＝Y₁-X₁，去除序列Z₁的首项和尾项，以生成单调振幅序列A＝{a(k)}，并生成序列Z₂＝Y₂-X₂，去除序列Z₂的首项和尾项，以生成单调周期序列C＝{c(k)}，其中，k＝1,2,3,…,n-1，X₁＝{0,p(i)}，Y₁＝{p(i),0}，X₁＝{0,t(i)}，Y₁＝{t(i),0}；

S3：将序列A和序列C组成向量序列O＝{(a(k),c(k))}，其中k＝1,2,3,…,n-1，振幅序列为A＝{a(k)}，周期序列为C＝{c(k)}，k＝1,2,3,…,n-1；

S4：从序列A中提取n-m+1个m维向量A(u)＝{a(u),a(u+1),a(u+2),…,a(u+m-1)}，从C中提取n-m+1个m维向量C(u)＝{c(u),c(u+1),c(u+2),…,c(u+m-1)}，并由A(u)和C(u)组成向量序列Q(u)＝{(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))}；

S5：计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离；

S6：计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]；

S7：设相似容限为R，统计出D[Q_m(u),Q_m(v)]<R的个数num{D[Q_m(u),Q_m(v)]<R}，计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)；

S8：计算所有向量B_u^m(R)的平均值B^m(R)；

S9：将维数由m加到m+1，重复上述步骤S4～S8，得到B^m+1(R)；

S10：根据所述B^m(R)和B^m+1(R)计算样本熵。

可选地，步骤S5中，通过下式计算计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离，

J_D((a(u),c(u)),(a(v),c(v)))＝1-J((a(u),c(u)),(a(v),c(v)))

其中，J_D((a(u),c(u)),(a(v),c(v)))为向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离，J((a(u),c(u)),(a(v),c(v)))为向量元素(a(u),c(u))和(a(v),c(v))的Jaccard相似系数，所述Jaccard相似系数的计算公式为ΔB为容纳(a(u),c(u))的最小矩形B(u)与容纳(a(v),c(v))的最小矩形B(v)的重合面积。

可选地，步骤S6中，通过下式计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]，

D[Q_m(u),Q_m(v)]＝max(J_D((a(u+x),c(u+x)),(a(v+x),c(v+x))))

其中，0≤x≤m-1，1≤u≤n-m+1，1≤v≤n-m+1。

可选地，步骤S7中，通过下式计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)，

$B_u^m\left(R\right)=\frac{1}{n-m}num\{D\&lsqb;Q_m\left(u\right),Q_m\left(v\right)\&rsqb;\&lt;R\}$

其中，1≤v≤n-m+1且v≠u。

可选地，步骤S10中，根据所述B^m(R)和B^m+1(R)通过下式计算样本熵，

D2SEn(m,R,n)＝-ln[B^m+1(R)/B^m(R)]

其中，D2SEn(m,R,n)为样本熵。

依据本发明的另一个方面，提供了一种基于二维空间的时间序列样本熵的计算系统，所述系统包括：

数据提取单元，用于提取时间序列的极值点和极值点的时间位置，以生成极值序列P＝{p(i)}及时间位置序列T＝{t(i)}，其中，i＝1,2,3,…,n；

序列生成单元，用于生成序列Z₁＝Y₁-X₁，去除序列Z₁的首项和尾项，以生成单调振幅序列A＝{a(k)}，并生成序列Z₂＝Y₂-X₂，去除序列Z₂的首项和尾项，以生成单调周期序列C＝{c(k)}，其中，k＝1,2,3,…,n-1，X₁＝{0,p(i)}，Y₁＝{p(i),0}，X₁＝{0,t(i)}，Y₁＝{t(i),0}；

序列组成单元，用于将序列A和序列C组成向量序列O＝{(a(k),c(k))}，其中k＝1,2,3,…,n-1，振幅序列为A＝{a(k)}，周期序列为C＝{c(k)}，k＝1,2,3,…,n-1；

向量提取单元，用于从序列A中提取n-m+1个m维向量A(u)＝{a(u),a(u+1),a(u+2),…,a(u+m-1)}，从C中提取n-m+1个m维向量C(u)＝{c(u),c(u+1),c(u+2),…,c(u+m-1)}，并由A(u)和C(u)组成向量序列Q(u)＝{(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))}；

元素计算单元，用于计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离；

向量计算单元，用于计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]；

比值计算单元，用于设相似容限为R，统计出D[Q_m(u),Q_m(v)]<R的个数num{D[Q_m(u),Q_m(v)]<R}，计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)；

均值计算单元，用于计算所有向量B_u^m(R)的平均值B^m(R)；

维数调整单元，用于将维数由m加到m+1，得到B^m+1(R)；

样本熵计算单元，用于根据所述B^m(R)和B^m+1(R)计算样本熵。

可选地，所述元素计算单元通过下式计算计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离，

J_D((a(u),c(u)),(a(v),c(v)))＝1-J((a(u),c(u)),(a(v),c(v)))

其中，J_D((a(u),c(u)),(a(v),c(v)))为向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离，J((a(u),c(u)),(a(v),c(v)))为向量元素(a(u),c(u))和(a(v),c(v))的Jaccard相似系数，所述Jaccard相似系数的计算公式为ΔB为容纳(a(u),c(u))的最小矩形B(u)与容纳(a(v),c(v))的最小矩形B(v)的重合面积。

可选地，所述向量计算单元通过下式计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]，

D[Q_m(u),Q_m(v)]＝max(J_D((a(u+x),c(u+x)),(a(v+x),c(v+x))))

其中，0≤x≤m-1，1≤u≤n-m+1，1≤v≤n-m+1。

可选地，所述比值计算单元通过下式计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)，

$B^m\left(R\right)=\frac{1}{n-m}num\{D\&lsqb;Q_m\left(u\right),Q_m\left(v\right)\&rsqb;\&lt;R\}$

其中，1≤v≤n-m+1且v≠u。

可选地，所述样本熵计算单元根据所述B^m(R)和B^m+1(R)通过下式计算样本熵，

D2SEn(m,R,n)＝-ln[B^m+1(R)/B^m(R)]

其中，D2SEn(m,R,n)为样本熵。

本发明通过各步骤之间的配合，实现了样本熵的计算，能够应用于脑电复杂度计算，也可应用于其他存在局部极值点的振动序列或波形的复杂度计算。振动熵可作为对波形进行模式识别的特征指标，对不同复杂度波形进行分类。在信号处理时，信号中混有白噪声的振幅和周期是在一定范围内呈混沌的正态分布的，而信号的波形较为规律，所以本发明也可用于对信号中噪声的识别和剔除。

附图说明

图1是本发明一种实施方式的基于二维空间的时间序列样本熵的计算方法的流程图；

图2是样本熵与相似容限R之间的关系图(其中，星号表示符合S1<S2<S3<S4<S5的关系)；

图3是在相似容限取一定值时，仿真波形的样本熵示意图；

图4是本发明一种实施方式的基于二维空间的时间序列样本熵的计算系统的结构框图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

图1是本发明一种实施方式的基于二维空间的时间序列样本熵的计算方法的流程图；参照图1，所述方法包括：

S1：提取时间序列的极值点和极值点的时间位置，以生成极值序列P＝{p(i)}及时间位置序列T＝{t(i)}，其中，i＝1,2,3,…,n；

S2：生成序列Z₁＝Y₁-X₁，去除序列Z₁的首项和尾项，以生成单调振幅序列A＝{a(k)}，并生成序列Z₂＝Y₂-X₂，去除序列Z₂的首项和尾项，以生成单调周期序列C＝{c(k)}，其中，k＝1,2,3,…,n-1，X₁＝{0,p(i)}，Y₁＝{p(i),0}，X₁＝{0,t(i)}，Y₁＝{t(i),0}；

S3：将序列A和序列C组成向量序列O＝{(a(k),c(k))}，其中k＝1,2,3,…,n-1，振幅序列为A＝{a(k)}，周期序列为C＝{c(k)}，k＝1,2,3,…,n-1；

S4：从序列A中提取n-m+1个m维向量A(u)＝{a(u),a(u+1),a(u+2),…,a(u+m-1)}，从C中提取n-m+1个m维向量C(u)＝{c(u),c(u+1),c(u+2),…,c(u+m-1)}，并由A(u)和C(u)组成向量序列Q(u)＝{(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))}；

需要说明的是，A(u)和C(u)为某段振动序列中的两个表示振动特征的向量。

S5：计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离；

可理解的是，在振幅-周期的欧式空间中，(a(u),c(u))与(a(v),c(v))的Jaccard相似系数为：容纳(a(u),c(u))的最小矩形B(u)与容纳(a(v),c(v))的最小矩形B(v)的重合面积ΔB，与此时B(u)和B(v)共占空间面积的比值为

$J\left(\right(a\left(u\right),c\left(u\right)),(a\left(v\right),c\left(v\right)\left)\right)=\frac{\mathrm{\&Delta;}B}{a\left(u\right)\&times;c\left(u\right)+a\left(v\right)\&times;c\left(v\right)-\mathrm{\&Delta;}B}$

其中ΔB＝min(a(u),a(v))×min(c(u),c(v))，是矩形B(u)和B(v)重合面积的最大值。因此，向量(a(u),c(u))与(a(v),c(v))的Jaccard距离计算公式为：

J_D((a(u),c(u)),(a(v),c(v)))＝1-J((a(u),c(u)),(a(v),c(v)))

S6：计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]；

其中，向量序列Q(u)与Q(v)的距离表示为：Q(u)与Q(v)序列对应向量元素Jaccard距离的最大值。

D[Q_m(u),Q_m(v)]＝max(J_D((a(u+x),c(u+x)),(a(v+x),c(v+x))))

其中，0≤x≤m-1，1≤u≤n-m+1，1≤v≤n-m+1。

S7：设相似容限为R，统计出D[Q_m(u),Q_m(v)]<R的个数num{D[Q_m(u),Q_m(v)]<R}，计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)；

在具体实现中，通过下式计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)，

$B^m\left(R\right)=\frac{1}{n-m}num\{D\&lsqb;Q_m\left(u\right),Q_m\left(v\right)\&rsqb;\&lt;R\}$

其中，1≤v≤n-m+1且v≠u。

S8：计算所有向量B_u^m(R)的平均值B^m(R)；

在具体实现中，通过下式计算所有向量B_u^m(R)的平均值B^m(R)，

$B^m\left(R\right)=\frac{1}{n-m+1}\underset{u=1}{\overset{n-m+1}{\&Sigma;}}{B}_u^m\left(R\right).$

S9：将维数由m加到m+1，重复上述步骤S4～S8，得到B^m+1(R)；

其中，

S10：根据所述B^m(R)和B^m+1(R)计算样本熵。

在具体实现中，根据所述B^m(R)和B^m+1(R)通过下式计算样本熵，

D2SEn(m,R,n)＝-ln[B^m+1(R)/B^m(R)]

其中，D2SEn(m,R,n)为样本熵。

可理解的是，为了便于对上述实施方式的方法进行验证，本实施方式中，在步骤S1之前，可依据波形特点由正弦波构造仿真波形，子波函数为：

$s_i\left(t\right)=A_i\&CenterDot;cos\left(\frac{2\mathrm{\&pi;}}{T_i}t\right)$

其中，当i＝1,3,5,…,2n-1时，t∈[0,π)，当i＝2,4,6,…,2n时，t∈[π,2π)，Ti和Ai分别表示第i个正弦子波的周期和振幅。

再设时间序列S的子波数为n，根据条件随机生成随机序列A＝{A₁,A₂,A₃,…,A_n}和周期序列T＝{T₁,T₂,T₃,…,T_n}，由A和T组成振动向量序列O＝{(A₁,T₁),(A₂,T₂),…,(A_n,T_n)}，然后由O中的各元素生成余弦子波序列，将序列中的子波收尾相连构成仿真时间序列S。

然后，构造五种波形S1、S2、S3、S4、S5进行实验。波形的复杂度关系为S1<S2<S3<S4<S5。

计算上述五种波形的样本熵，振动熵的计算方法则按照上述步骤S1～S10进行计算，当然，在进行计算前，先对上述模拟信号进行滤波，将其分解为若干子频段信号。

由图2可知，D2SEn的值与相似容限R的取值有关，D2SEn随相似容限R的增大而减小，仅当相似容限的取值在一定范围内时，仿真波形D2SEn值的大小关系才与预设条件相同。

图3指示在相似容限取一定值时仿真波形的二维样本熵示意图。分别对R＝0.05、0.1、0.15时的D2SEn值进行方差分析。结果表明，在上述三种条件下，D2SEn值在五种波形间的组效应显著(R＝0.05：F＝3245.084,P<0.00001；R＝0.1：F＝3245.084,P<0.00001；R＝0.15：F＝3923.475,P<0.00001)。组间两两比较时，每两组波形间的D2SEn值均存在显著差异(P<0.00001)。由图2可知，当相似容限R的取值在一定范围内时，D2SEn值可以有效反映波形的复杂度特征。

本实施方式的方法可应用于脑电复杂度计算，也可应用于其他存在局部极值点的振动序列或波形的复杂度计算。振动熵可作为对波形进行模式识别的特征指标，对不同复杂度波形进行分类。在信号处理时，信号中混有白噪声的振幅和周期是在一定范围内呈混沌的正态分布的，而信号的波形较为规律，所以本实施方式的方法也可用于对信号中噪声的识别和剔除。

图4是本发明一种实施方式的基于二维空间的时间序列样本熵的计算系统的结构框图；参照图4，所述系统包括：

数据提取单元401，用于提取时间序列的极值点和极值点的时间位置，以生成极值序列P＝{p(i)}及时间位置序列T＝{t(i)}，其中，i＝1,2,3,…,n；

序列生成单元402，用于生成序列Z₁＝Y₁-X₁，去除序列Z₁的首项和尾项，以生成单调振幅序列A＝{a(k)}，并生成序列Z₂＝Y₂-X₂，去除序列Z₂的首项和尾项，以生成单调周期序列C＝{c(k)}，其中，k＝1,2,3,…,n-1，X₁＝{0,p(i)}，Y₁＝{p(i),0}，X₁＝{0,t(i)}，Y₁＝{t(i),0}；

序列组成单元403，用于将序列A和序列C组成向量序列O＝{(a(k),c(k))}，其中k＝1,2,3,…,n-1，振幅序列为A＝{a(k)}，周期序列为C＝{c(k)}，k＝1,2,3,…,n-1；

向量提取单元404，用于从序列A中提取n-m+1个m维向量A(u)＝{a(u),a(u+1),a(u+2),…,a(u+m-1)}，从C中提取n-m+1个m维向量C(u)＝{c(u),c(u+1),c(u+2),…,c(u+m-1)}，并由A(u)和C(u)组成向量序列Q(u)＝{(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))}；

元素计算单元405，用于计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离；

向量计算单元406，用于计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]；

比值计算单元407，用于设相似容限为R，统计出D[Q_m(u),Q_m(v)]<R的个数num{D[Q_m(u),Q_m(v)]<R}，计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)；

均值计算单元408，用于计算所有向量B_u^m(R)的平均值B^m(R)；

维数调整单元409，用于将维数由m加到m+1，得到B^m+1(R)；

样本熵计算单元410，用于根据所述B^m(R)和B^m+1(R)计算样本熵。

可选地，所述元素计算单元通过下式计算计算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离，

J_D((a(u),c(u)),(a(v),c(v)))＝1-J((a(u),c(u)),(a(v),c(v)))

其中，J_D((a(u),c(u)),(a(v),c(v)))为向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距离，J((a(u),c(u)),(a(v),c(v)))为向量元素(a(u),c(u))和(a(v),c(v))的Jaccard相似系数，所述Jaccard相似系数的计算公式为ΔB为容纳(a(u),c(u))的最小矩形B(u)与容纳(a(v),c(v))的最小矩形B(v)的重合面积。

可选地，所述向量计算单元通过下式计算向量序列Q(u)与Q(v)的距离D[Q_m(u),Q_m(v)]，

D[Q_m(u),Q_m(v)]＝max(J_D((a(u+x),c(u+x)),(a(v+x),c(v+x))))

其中，0≤x≤m-1，1≤u≤n-m+1，1≤v≤n-m+1。

可选地，所述比值计算单元通过下式计算num{D[Q_m(u),Q_m(v)]<R}与所有距离总数的比值B_u^m(R)，

$B_u^m\left(R\right)=\frac{1}{n-m}num\{D\&lsqb;Q_m\left(u\right),Q_m\left(v\right)\&rsqb;\&lt;R\}$

其中，1≤v≤n-m+1且v≠u。

可选地，所述样本熵计算单元根据所述B^m(R)和B^m+1(R)通过下式计算样本熵，

D2SEn(m,R,n)＝-ln[B^m+1(R)/B^m(R)]

其中，D2SEn(m,R,n)为样本熵。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。