两种滑动排列熵序列突变检测的对比分析

罗文相，万 丽，2，赖斯敏

(1.广州大学数学与信息科学学院，广东 广州 510006; 2.广州大学数学与交叉科学广东普通高校重点实验室，广东 广州 510006)

随着非线性科学的发展，时间序列突变检测越来越得到人们的关注。传统的序列突变检测方法，检测结果严重依赖于时间序列的长度[1]。因此，许多学者将熵的理论引入到时间序列突变检测领域。熵是系统复杂度的一种度量，可用于时间序列动力学结构突变的检测。在信息论建立后，关于熵的理论和应用得到了迅速发展[2]。Pincus等[3-4]基于边缘概率分布统计衡量时间序列复杂度提出的近似熵,被广泛地应用于生理序列分析、矿化识别及径流突变分析等领域[5-7]; 针对近似熵存在自身匹配特性的缺点，Richman等[8]对近似熵作了改进，提出了样本熵。进一步，Bandt等[9]提出了排列熵(permutation entropy，PE)，与近似熵、样本熵及分形维数、Lyapunov指数等复杂度参数相比，具有算法简单，抗噪能力强等特点，能有效地放大序列的微小变化，广泛应用于信号突变检测，并在机械故障诊断、医学和气候等领域取得良好的应用效果[10-15]。

目前，对排列熵的应用大多只是单纯地计算出时间序列的单个熵值，根据熵值的大小来判断序列的复杂度。本文将排列熵与滑动窗口和滑动移除数据技术相结合，比较了滑动排列熵(moving permutation entropy，M-PE) 和滑动移除排列熵 (moving cut data-permutation entropy，MC-PE ) 方法在非线性时间序列动力学结构突变检测的性能，为排列熵的实际应用提供更好的参考方法[16-17]。

1 时间序列突变检测方法——排列熵

1.1 PE方法

PE算法的基本原理如下:

设一长度为N的时间序列X= {x(n)，n=1，2,…，N}，对时间序列X进行相空间重构，得到N- (m- 1)τ行，m列的矩阵

Xk=[x(i)，x(i+τ),…，x(i+ (m-1)τ)] ，i=1，2,…，N-(m-1)τ。

(1)

(1)式中：m是嵌入维，τ是延迟时间。Xk中的每一行为一个重构分量，共有N-(m-1)τ个。将第i个重构分量 (x(i)，x(i+τ),…，x(i+ (m- 1)τ)) 按照升序重新排列，得到

x(i+(j1-1)τ) ≤x(i+ (j2-1)τ) ≤…≤x(i+ (jm- 1)τ)。

(2)

(2)式中：j1、j2、…、jm表示元素所在重构分量中列的索引，若元素相等，则按照元素索引的大小进行排列。因此，矩阵Xk中的每个重构分量重新按照升序排列，可得到一个N-(m-1)τ行m列的符号序列矩阵A。记A(i)=(j1,j2,…，jm)，为(1，2,…，m)的某一种排列，显然m个元素最多有m!种不同排列。以时间序列{5，3，8，6，4}为例，当m=3，τ= 1时，得到的第一个重构分量为(5，3，8)，由于xt+1

对所有A(i)的排列情况进行统计，记录A(i)中某种排列A(r)的个数为N(r)，其中r≤m!，0

(3)

则排列熵(E)的定义为:

E(m)=-∑P(r)log(P(r))。

(4)

(4)式中：log表示为以2为底的对数。当P(r)=1/(m!) 时，E(m)取得最大值log(m!)。在实际应用中，通常会对E(m)进行归一化处理，即:

0 ≤E(m)=E(m)/log(m!) ≤ 1。

(5)

E(m)值的大小与时间序列X的随机程度表现一致，E(m)值越大，说明时间序列X的随机程度越高，反之则X越规则。

1.2 M-PE方法

M-PE方法是一种将排列熵与滑动窗口技术相结合的一种突变检测方法，步骤是对一时间序列，以长度为W的滑动窗口选取时间序列的数据构成子序列，采用PE方法计算该子序列的排列熵值，保持滑动窗口W不变，以步长为S逐步滑动选取新的子序列，同样采用PE方法计算新子序列的排列熵值，直到时间序列结束，得到一个随步长逐步变化的排列熵值序列。根据排列熵值序列的波动变化，判断序列的突变点或突变区间。

1.3 MC-PE方法

MC-PE方法是一种将排列熵与滑动移除数据技术相结合的一种突变检测方法,具体的步骤如下:

1) 选择滑动移除数据的窗口长度M；

2) 从时间序列{x(n)}的第t(t=1，2,…，N-M+ 1，N为时间序列{x(n)}的长度)个数据开始连续移除M个数据，再将剩余的N-M个数据连在一起，得到一个长度为N-M的新子序列；

3) 利用排列熵方法计算新子序列的排列熵值；

4) 保持数据移除的窗口长度不变，取滑动移除步长为M，逐步移动窗口，重复2～3步操作，可得到一个长度为int(N/M)(int表示取整)的排列熵序列；

5) 基于相同动力学性质的数据复杂度差异不大，而不同动力学性质的数据复杂度不相同这一特点，结合步骤4得到的排列熵序列确定突变点或突变区间。

对于同一个具有相关性时间序列，任意移除动力学性质相同的数据，对其E(m)值的影响几乎可以忽略不计。因此，可以通过步骤1～5观察E(m)值的变化来检测时间序列的动力学结构突变。

2 M-PE和MC-PE方法的比较

2.1 非线性理想时间序列的构造

构造由Logistic映射产生的非线性理想时间序列IS，序列总长度为2 000，其方程如下:

xn+1=λxn(1-xn)，xn∈[0,1] 。

(6)

图1 非线性理想时间序列ISFig. 1 Nonlinear ideal time series IS

(6)式中：初始值为x0= 0.7，λ∈[0,4] ，表示虫口模型的生长率，IS的前1 000个数据的参数为λ= 3.6，后1 000个数据的参数变为λ= 3.7，两者均处于混沌状态，其演化曲线由图1给出。从IS的构造可知，在n=1 001时，系统的演化由较低的生长率3.6变为较高的生长率3.7。排列熵算法的参数均取m=3,τ=2,得出序列IS前1 000个数据的排列熵值为0.575 2；后1 000个数据的排列熵值为0.661 2。IS后1 000个数据的排列熵值比前1 000个数据大，由排列熵的物理意义可知,在n=1 001后系统的复杂度增大，可预测性变小。

2.2 M-PE的非线性时间序列状态识别

M-PE方法中滑动窗口W需要选取适当的大小，如果滑动窗口过小，则因所含数据量小，不能很好地反应原始序列的特征；如果滑动窗口过大，则得到的滑动排列熵值个数少，忽略了原始序列的部分细致特点。经过多组数值试验，发现选取大小为原始序列长度的5%～10%的滑动窗口能取得比较好的效果。选取滑动窗口W=150,对理想时间序列IS进行滑动排列熵的计算，滑动步长分别为S=10，20，30和50的检测结果由图2给出。

图2 非线性时间序列IS的M-PE检测结果Fig. 2 M-PE detection results of nonlinear time series IS

由图2可见，对不同的滑动步长，M-PE方法得到的E(m)值变化趋势基本一致，在n=1 000前后，E(m)值呈现出两种不同的稳定状态。前半部分的E(m)值处于较低水平，而后半部分的E(m)值处于较高水平，由排列熵物理意义，可知前半部分数据的复杂度比后半部分的复杂度低，与理想时间序列IS的构造基本一致，说明M-PE方法具有识别序列动力学结构突变的能力。从图2可以清楚地看到，在前后不同水平的E(m)值中间存在一个长度与滑动窗口相近的上升区域，这是因为随着滑动窗口的移动，滑动窗口中所包含的数据由完全来自IS前1 000的数据，逐渐变为前1 000数据和后1 000数据的混合，最后只含有来自IS后1 000的数据。区间(1 000，1 150)是两种数据混合的阶段，在这个阶段，由于后1 000数据的个数逐渐增加，滑动窗口数据的复杂度也随之增加，直到滑动窗口数据的复杂度接近于后1 000数据的复杂度，因此形成了与滑动窗口大小相近的上升突变区间。从上升区间的分析中，可以判断突变点的位置大致在上升区间的开始处，即约在n=1 001处，说明M-PE方法能检测时间序列的动力学结构突变。

2.3 MC-PE的非线性时间序列状态识别

作为对比，MC-PE方法中滑动移除步长(即滑动移除窗口)同样分别选择M= 10，20，30和50，其检测结果如图3所示。

图3 非线性时间序列IS的MC-PE检测结果Fig.3 MC-PE detection results of nonlinear time series IS

4种滑动移除步长得到的MC-PE检测结果几乎是一致的，以n=1 000为界，移除相应数据后得到的排列熵值分为两种不同的演化阶段，在n>1 000的E(m)值明显小于n≤1 000时的情形，而前后两种阶段的E(m)值各自处于相同的状态。根据排列熵的物理意义，排列熵值越小表示时间序列的复杂度越小，而移除数据后得到较小的E(m)值意味着所移除的数据的复杂度较大，即移除复杂度较大的数据使得剩余数据的复杂度变小，说明序列在n>1 000处的复杂度高于前半部分的复杂度，这与理想时间序列IS的构造完全一致。图3(a)～3(d)均显示出E(m)值在n=1 000前后存在明显差异，而且随着滑动移除窗口M的增大，移除前后不同动力学结构的数据后所得的E(m)值之间的差异逐渐增大，因此可以判断IS的突变点大约在n=1 001处，说明MC-PE方法能有效地识别时间序列的动力学结构突变。

2.4 M-PE和MC-PE结果对比分析

M-PE和MC-PE方法均能识别时间序列的动力学结构突变，但M-PE方法的结果存在一个大小与滑动窗口相近的上升突变区间，当滑动窗口增大时，上升区间也随之增大，因此比较难精确地判定上升突变区间的开始位置，即M-PE方法得出的突变点位置误差比较大。另外，滑动窗口W的选取与时间序列的长度有关,选取滑动窗口W=30和50，滑动步长S=1对IS进行M-PE检测，结果显示，在n=1 000附近的突变区域分界不明确，而且前半部分和后半部分各自的波动比较大(图略)，不能有效地判断突变点的位置。因此，不能为了提高M-PE突变检测结果的精确度而刻意减少滑动窗口W的大小。

相对于M-PE方法依赖于滑动窗口的大小，导致突变点位置误差较大这一缺点，MC-PE方法的突变检测结果就更为精确。从图3的结果可以看出，前半部分和后半部分的E(m)值之间的突变区域与滑动移除窗口M(即滑动移除步长)相近，换言之，MC-PE方法得出的突变点位置可以精确到滑动移除窗口M范围内，这大大提高了突变点的精确度。另外，MC-PE方法几乎不依赖于滑动移除窗口的大小，对于选取较小的滑动移除窗口，MC-PE方法依然能识别出突变点的位置。图4给出了滑动移除窗口M=1和3时，MC-PE方法对时间序列IS的突变检测结果。

图4 IS的MC-PE检测结果Fig.4 MC-PE detection results of IS

如图4(a)所示,可以十分清晰地看到E(m)值以n=1 000为界，前后分别处于两种不同的状态，由于M=1，因此可以精确地判断IS的突变点为n=1 001。在实际中，一般不会选取M=1，而是选取相对大一点的滑动移除窗口，因为每次只移除一个数据，计算成本太高，而且只移除一个数据，对实测数据的影响不是十分明显。图4(b)的结果表明，滑动移除窗口M=3比M=1更能反映序列IS后1 000数据的复杂度高于前1 000数据的复杂度，因为图4(b)后半部分的E(m)值绝大部分都处于较低水平。综上所述，M-PE方法对时间序列突变具有一定的识别能力，但该方法依赖于滑动窗口的选取，阻碍了其在实际中的应用。 MC-PE方法的检测结果几乎不依赖于滑动移除窗口的大小，能够更为精准地得出时间序列的突变位置，明显优于M-PE方法。

3 结语

基于衡量时间序列复杂度的参数——排列熵，结合滑动窗口和滑动移除数据技术，比较了M-PE和MC-PE方法在非线性时间序列动力学结构突变检测的性能。分析表明：M-PE方法对时间序列突变具有一定的识别能力，但检测结果存在一个大小与滑动窗口相近的上升突变区间，得出的突变点位置误差比较大，且该方法依赖于滑动窗口的选取，阻碍了在实际中的应用。与M-PE相比，MC-PE方法的检测结果几乎不依赖于滑动移除窗口的大小，能够更为精准地得出时间序列的突变位置，即便是选择非常小的滑动移除窗口 (M=1)，MC-PE方法也能有效地识别出序列IS的突变点，明显优于M-PE方法，具有良好的应用前景。