基于樣本熵的時(shí)間序列非線性檢測方法研究

李聰改，曹 銳，武 政，相 潔，2

（1.太原理工大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，山西 太原030024；2.北京工業(yè)大學(xué) 國際WIC研究院，北京100000）

0 引 言

對于生理信號等觀測得到的時(shí)間序列，目前常用于檢驗(yàn)其非線性特性的方法主要是替代數(shù)據(jù)法。替代數(shù)據(jù)法最早是在1992年由Theiler等人提出，之后很多學(xué)者對該方法提出了改進(jìn)［3－9］。用于驗(yàn)證原數(shù)據(jù)和替代數(shù)據(jù)之間差別的特 征 量 有 很 多，如 關(guān) 聯(lián) 維D2［10－12］、近 似 熵 （approximate entropy，ApEn）［13－16］等。然而，這些指標(biāo)都存在一些不足，如使用關(guān)聯(lián)維時(shí)，只有數(shù)據(jù)量很大時(shí)才能得到可靠結(jié)論［17］。近似熵的結(jié)論也缺乏相對一致性［18］，且計(jì)算時(shí)間較長。樣本熵 （sample entropy，SampEn）是近似熵算法的改進(jìn)，比近似熵具有更好的穩(wěn)定性。樣本熵是否可以作為特征量來判定時(shí)間序列的非線性特性？本文針對這一問題進(jìn)行了相關(guān)研究，文中將樣本熵作為替代數(shù)據(jù)法中的特征量，在3種仿真時(shí)間序列上驗(yàn)證樣本熵是否可以作為特征量來識別時(shí)間序列的非線性特性。本文同時(shí)也記錄了近似熵作為特征量的實(shí)驗(yàn)結(jié)果，并對近似熵和樣本熵的結(jié)果進(jìn)行了對比。

1 基于樣本熵的替代數(shù)據(jù)法

1.1 替代數(shù)據(jù)分析方法

替代數(shù)據(jù)法中，首先提出零假設(shè) （假設(shè)原始待測數(shù)據(jù)是線性的），其次通過振幅調(diào)節(jié)傅里葉變換 （AAFT）等方法產(chǎn)生替代數(shù)據(jù)。替代數(shù)據(jù)是隨機(jī)排列的線性數(shù)據(jù)，它與原始數(shù)據(jù)具有相同的功率譜，但是，原始數(shù)據(jù)中的非線性成分已經(jīng)被完全破壞掉了。產(chǎn)生替代數(shù)據(jù)之后，分別計(jì)算原始數(shù)據(jù)與替代數(shù)據(jù)的特征量，并采用如式 （1）所示的統(tǒng)計(jì)檢驗(yàn)方法來比較原始數(shù)據(jù)和替代數(shù)據(jù)之間的差異，如果原始數(shù)據(jù)與替代數(shù)據(jù)特征量差別顯著，則拒絕原假設(shè)，即原始時(shí)間序列為非線性；反之為線性

式中：Ao——原始數(shù)據(jù)的某特征量的值，As——多組替代數(shù)據(jù)的該特征量的平均值，SD（As）——多組替代數(shù)據(jù)的該特征量的標(biāo)準(zhǔn)差。

1.2 樣本熵

樣本熵SampEn是由Richman和Moorman提出的一種新的時(shí)間序列復(fù)雜性測度方法［19］，是ApEn的一種改進(jìn)算法，其計(jì)算方法如下：

步驟1 將序列x（1），x（2），...，x（N）按順序組成m維矢量，即

步驟2 定 義Xm（i）與Xm（i）間 的 距 離d［Xm（i），Xm（i）］為兩者對應(yīng)元素中差值最大的一個(gè)，即

式中：1≤k≤m－1；1≤i，j≤N－m＋1，i≠j。

步驟3 給定的閾值r（r＞0），對每一個(gè)1≤i≤Nm 值統(tǒng)計(jì)d［Xm（i），Xm（i）］＜r 的數(shù)目 （稱之為模板匹配數(shù)）與總的距離數(shù)目N－m－1的比值，記作Bmi（r）

其中1≤j≤N－m，j≠i。求其對i所有的平均值

對于m＋1點(diǎn)矢量，同樣有

其中1≤j≤N－m，j≠i。求其對所有i的平均值

步驟4 此序列的樣本熵為

但是在實(shí)際計(jì)算中N 不可能為∞，當(dāng)N 取有限值時(shí)，估計(jì)

樣本熵的算法思想是：當(dāng)維數(shù)為m 時(shí)，時(shí)間序列的自相似概率用B 表示；當(dāng)維數(shù)為m＋1 時(shí)，時(shí)間序列的自相似概率用A 表示，從而可以得出公式：CP＝A／B?？梢钥闯觯瑯颖眷氐挠?jì)算是以－ln （CP）為模型。近似熵的算法與樣本熵類似，但是近似熵的計(jì)算過程中對數(shù)據(jù)自身進(jìn)行了匹配，因而在計(jì)算的過程中肯定會存在偏差。另外近似熵的值與相似容限r(nóng)的取值有很大關(guān)系，導(dǎo)致對于不同的參數(shù)計(jì)算結(jié)果缺乏相對一致性。因此本文提出將近似熵的改進(jìn)算法—樣本熵作為特征量引入到替代數(shù)據(jù)法中來檢驗(yàn)時(shí)間序列的非線性特征。

2 仿真數(shù)據(jù)

本文實(shí)驗(yàn)所采用的數(shù)據(jù)為線性AR 模型、Lorenz方程、Logistic方程所產(chǎn)生的時(shí)間序列。其中線性AR 模型產(chǎn)生的時(shí)間序列確定是線性的；Lorenz方程產(chǎn)生的是典型的連續(xù)的非線性時(shí)間序列；而Logistic方程產(chǎn)生的是典型的離散的非線性時(shí)間序列。實(shí)驗(yàn)將驗(yàn)證本文提出的基于樣本熵的非線性檢測方法對這3種時(shí)間序列的檢測結(jié)果是否與時(shí)間序列原特性一致。

為了驗(yàn)證SampEn作為特征量時(shí)，是否如關(guān)聯(lián)維D2一樣會受到時(shí)間序列長度的影響，實(shí)驗(yàn)中采用8種不同的時(shí)間長度的序列。同時(shí)，為了避免計(jì)算結(jié)果的偶然性，仿真實(shí)驗(yàn)為每個(gè)時(shí)間長度均產(chǎn)生了20組原始數(shù)據(jù)。另外，實(shí)驗(yàn)中記錄了SampEn與ApEn作為特征量進(jìn)行非線性檢測時(shí)的運(yùn)算時(shí)間，并進(jìn)行了對比。

2.1 非線性數(shù)據(jù)

本文中采用兩種方程來產(chǎn)生確定性混沌的時(shí)間序列：利用Lorenz方程產(chǎn)生連續(xù)的非線性時(shí)間序列，用Logistic方程產(chǎn)生離散的非線性時(shí)間序列。

Lorenz方程公式如下所示

當(dāng)參數(shù)σ＝10；r＝28；b＝8／3時(shí)，方程產(chǎn)生的數(shù)據(jù)為非線性。

產(chǎn)生非線性數(shù)據(jù)的Logistic離散映射方程如下所示

當(dāng)μ＝4.0，初始值x （1）＝0.1的時(shí)候，產(chǎn)生的數(shù)據(jù)為非線性。

2.2 線性數(shù)據(jù)

本文采用一階自回歸AR 模型來產(chǎn)生線性數(shù)據(jù)，其計(jì)算方法如下所示

式中：e（t）——白噪聲，服從均值為0，方差σ2的分布。論文中，取系數(shù)k＝0.2，取初始值y（1）＝0，e（t）的均值為0，方差為1。

為了避免數(shù)據(jù)的隨機(jī)性，實(shí)驗(yàn)中首先對上述3種待測數(shù)據(jù)分別產(chǎn)生長度為5000 的數(shù)據(jù)，分別隨機(jī)選擇長度為100，200，300，500，800，1000，2000，3000 的數(shù)據(jù) 各20次作為原始數(shù)據(jù)并計(jì)算特征量，其次對每個(gè)長度下每組原始數(shù)據(jù)分別產(chǎn)生20組替代數(shù)據(jù)，并求得20組替代數(shù)據(jù)的特征量，最后計(jì)算得到原始數(shù)據(jù)和替代數(shù)據(jù)特征量的差別顯著度。ApEn 和SampEn 的計(jì)算過程中，涉及到參數(shù)m，r的選取。根據(jù)參考文獻(xiàn) ［19］，本文計(jì)算過程中，維數(shù)m＝2，相似容限r(nóng)＝0.15SD （SD 為時(shí)間序列的標(biāo)準(zhǔn)差）。

3 實(shí)驗(yàn)結(jié)果及分析

3.1 Lorenz

實(shí)驗(yàn)計(jì)算了8 種長度數(shù)據(jù)的特征量，圖1 （a）、圖1（b）僅給出數(shù)據(jù)長度N 分別為300和800時(shí)第1次產(chǎn)生的原始數(shù)據(jù)和對應(yīng)的20組替代數(shù)據(jù)的樣本熵值分布情況，其它長度結(jié)果類似。其中，橫坐標(biāo)sur表示第幾組替代數(shù)據(jù)，縱坐標(biāo)SampEn表示該數(shù)據(jù)的樣本熵值?！埃北硎驹紨?shù)據(jù)的樣本熵，“＊”表示替代數(shù)據(jù)的樣本熵。從圖1中可以看出，本次原始數(shù)據(jù)的熵值和替代數(shù)據(jù)的熵值存在差異。所有數(shù)據(jù)的差別顯著度結(jié)果見表1。

圖1 Lorenz時(shí)間序列樣本熵結(jié)果分布

表1 Lorenz方程產(chǎn)生數(shù)據(jù)的熵值比較（df＝19，ɑ＝95%）

從表1中可以看出，所有差別顯著度值均大于2.093（df＝19，ɑ＝0.95時(shí)的t檢驗(yàn)臨界值），意味著所有的差別顯著度均達(dá)到了統(tǒng)計(jì)檢驗(yàn)顯著性，無論近似熵還是樣本熵都可以判斷出原始時(shí)間序列中包含非線性成分 （和Lorenz時(shí)間序列已知特性相符）。此外，不同長度下非線性檢測的結(jié)果穩(wěn)定，表明非線性檢測結(jié)果不受數(shù)據(jù)長度的影響，并且樣本熵的計(jì)算用時(shí)比近似熵大幅度減少。

3.2 Logistic離散映射

對于Logitic方程產(chǎn)生的時(shí)間序列，同樣進(jìn)行了相關(guān)處理并給出了實(shí)驗(yàn)結(jié)果。圖2 （a）、圖2 （b）也僅給出數(shù)據(jù)長度N 分別為300和800時(shí)第1次產(chǎn)生的原始數(shù)據(jù)和對應(yīng)的20組替代數(shù)據(jù)的樣本熵值分布情況。所有數(shù)據(jù)的差別顯著度結(jié)果見表2。

圖2 Logistic時(shí)間序列樣本熵結(jié)果分布

表2 Logistic方程產(chǎn)生數(shù)據(jù)的熵值比較（df＝19，ɑ＝95%）

從表2中可以看出，Logistic數(shù)據(jù)與Lorenz數(shù)據(jù)的實(shí)驗(yàn)結(jié)果類似，所有的差別顯著度值均大于2.093 （df＝19，α＝0.95時(shí)的t檢驗(yàn)臨界值），表明樣本熵作為特征值同樣可以檢測出待測數(shù)據(jù)中包含了非線性的成分 （和Logistic時(shí)間序列已知特性相符），并且不受數(shù)據(jù)長度的影響，以及計(jì)算用時(shí)比近似熵短。

3.3 一階自回歸AR 模型

圖3 線性AR 時(shí)間序列樣本熵結(jié)果分布

同前兩個(gè)實(shí)驗(yàn)一樣，本次也計(jì)算了AR 模型產(chǎn)生的8種長度數(shù)據(jù)的特征量，圖3 （a）、圖3 （b）僅給出數(shù)據(jù)長度N 分別為300和800時(shí)第1次產(chǎn)生的原始數(shù)據(jù)和對應(yīng)的20組替代數(shù)據(jù)的樣本熵值分布情況，其它長度結(jié)果類似。但圖3 （a）、圖3 （b）中可以看出，本次20組替代數(shù)據(jù)的熵值分布在原始數(shù)據(jù)的熵值周圍，并沒有出現(xiàn)像圖1和圖2中熵值的分離情況，意味著這組原始數(shù)據(jù)與替代數(shù)據(jù)并沒有顯著的差異。所有數(shù)據(jù)的差別顯著度結(jié)果見表3。

表3 一階自回歸模型產(chǎn)生數(shù)據(jù)的兩種熵值比較（df＝19，ɑ＝95%）

表3中所有的差別顯著度值均小于2.093 （df＝19，α＝0.95時(shí)的t檢驗(yàn)臨界值），表明原始數(shù)據(jù)與替代數(shù)據(jù)近似熵與樣本熵的值的差異沒有達(dá)到統(tǒng)計(jì)檢驗(yàn)顯著性，即原始數(shù)據(jù)和替代數(shù)據(jù)不存在顯著差異，可以判斷該時(shí)間序列中不包含非線性成分，是線性的 （和AR 時(shí)間序列已知特性相符），并且該結(jié)果不受數(shù)據(jù)長度的影響，另外樣本熵計(jì)算用時(shí)比近似熵要短。

在本實(shí)驗(yàn)中，通過替代數(shù)據(jù)法對不同長度的時(shí)間序列進(jìn)行了非線性檢測。實(shí)驗(yàn)結(jié)果表明，近似熵和樣本熵這兩個(gè)指標(biāo)均可作為替代數(shù)據(jù)法的特征量，較好的驗(yàn)證時(shí)間序列中的非線性成分。但是在相同長度的時(shí)間序列上，樣本熵的用時(shí)要比近似熵短。

基于關(guān)聯(lián)維的時(shí)間序列非線性檢測對數(shù)據(jù)長度要求較高，本實(shí)驗(yàn)通過8種數(shù)據(jù)長度分析了樣本熵作為特征量檢驗(yàn)非線性特征時(shí)是否依賴于數(shù)據(jù)長度。結(jié)果表明樣本熵作為特征量時(shí)對數(shù)據(jù)長度要求不高，適合較短時(shí)間序列的非線性特性檢測。

4 結(jié)束語

論文在三組仿真數(shù)據(jù)上，驗(yàn)證了樣本熵是否可以作為特征量來驗(yàn)證時(shí)間序列中有無非線性成分。實(shí)驗(yàn)結(jié)果表明，對于不同長度的時(shí)間序列，樣本熵均可以作為替代數(shù)據(jù)法的特征量檢驗(yàn)出時(shí)間序列中是否含有非線性成分；并且采用樣本熵的替代數(shù)據(jù)法比采用近似熵的替代數(shù)據(jù)法計(jì)算效率得到明顯提高。所以，本研究得到結(jié)論：基于樣本熵的替代數(shù)據(jù)法是一種準(zhǔn)確、高效的非線性檢測方法。當(dāng)然本研究還需要進(jìn)一步在其它數(shù)據(jù)上進(jìn)行驗(yàn)證，來支持論文的結(jié)論。但是本文的結(jié)論對于時(shí)間序列的非線性檢測，以及非線性動(dòng)力學(xué)在時(shí)間序列分析中的應(yīng)用，特別是越來越多的生理信號的分析，有著較大的價(jià)值。

［1］Jouny C C，Bergey G K.Characterization of early partial seizure onset：Frequency，complexity and entropy ［J］.Clinical Neurophysiology，2012，123 （4）：658－669.

［2］Chen Z，Brown E N，Barbieri R.Characterizing nonlinear heartbeat dynamics within a point process framework ［J］.IEEE Transactions on Biomedical Engineering，2010，57（6）：1335－1347.

［3］Faes L，Zhao H，Chon K H，et al.Time－varying surrogate data to assess nonlinearity in nonstationary time series：Application to heart rate variability ［J］.IEEE Transactions on Biomedical Engineering，2009，56 （3）：685－695.

［4］Suzuki T，Ikeguchi T，Suzuki M.Algorithms for generating surrogate data for sparsely quantized time series ［J］.Physica D：Nonlinear Phenomena，2007，231 （2）：108－115.

［5］Keylock C J.A wavelet－based method for surrogate data generation ［J］.Physica D：Nonlinear Phenomena，2007，225（2）：219－228.

［6］Lucio JH，Valdes R，Rodriguez LR.Improvements to surrogate data methods for nonstationary times series ［J］.Physica Review E，2012，85 （5）：056202.

［7］Spasi c＇S.Surrogate data test for nonlinearity of the rat cerebellar electrocorticogram in the model of brain injury ［J］.Signal Processing，2010，90 （12）：3015－3025.

［8］Suzuki T，Ikeguchi T，Suzuki M.Algorithms for gene－rating surrogate data for sparsely quantized time series ［J］.Physica D：Nonlinear Phenomena，2007，231 （2）：108－115.

［9］Keylock C J.A wavelet－based method for surrogate data generation ［J］.Physica D：Nonlinear Phenomena，2007，225（2）：219－228.

［10］Spasic Sladana.Surrogate data test for nonlinearity of the rat cerebellar electrocorticogram in the model of brain injury ［J］.Signal Processing，2010，90 （12）：3015－3025.

［11］HUANG Xiaona，ZHANG Junying.The nonlinear analysis of epilepsy EEG signals［J］.Editorial Department of Hexi University，2010，26 （2）：49－52 （in Chinese）.［黃小娜，張軍英.癲癇腦電信號的非線性分析方法 ［J］.河西學(xué)院學(xué)報(bào)，2010，26 （2）：29－32.］

［12］WANG Liyuan.The pathological fetal heart rate signal analysis used surrogate data ［J］.Journal of Changchun Normal University，2007，26 （2）：49－52 （in Chinese）. ［王立媛.病態(tài)胎兒心率信號的替代數(shù)據(jù)分析 ［J］.長春師范學(xué)院學(xué)報(bào)（自然科學(xué)報(bào)），2007，26 （2）：49－52.］

［13］Sohn H，Kim I，Lee W，et al.Linear and non－linear EEG analysis of adolescents with attention－deficit／hyperactivity disorder during a cognitive task ［J］.Clinical Neurophysiology，2010，121 （11）：1863－1870.

［14］Ocak H.Automatic detection of epileptic seizures in EEG using discrete wavelet transform and approximate entropy［J］.Expert Systems with Applications，2009，36 （2）：2027－2036.

［15］Sohn H，Kim I，Lee W，et al.Linear and non－linear EEG analysis of adolescents with attention－deficit／hyperactivity disorder during a cognitive task ［J］.Clinical Neurophysiology，2010，121 （11）：1863－1870.

［16］Lodha N，Naik S K，Coombes S A，et al.Force control and degree of motor impairments in chronic stroke ［J］.Clinical Neurophysiology，2010，121 （11）：1952－1961.

［17］YU Qing.The analysis of correlation dimension calculation［J］.Journal of Tianjin University of Technology，2004，20（4）：60－62 （in Chinese）.［于青.關(guān)聯(lián)維數(shù)計(jì)算的分析研究［J］.天津理工學(xué)院學(xué)報(bào)，2004，20 （4）：60－62.］

［18］Cao Rui，Li Li，Chen Junjie.Comparative study of approximate entropy and sample entropy in EEG data analysis［J］.BioTechnology：An Indian Journal，2013，7（11）：493－498.

［19］Moorman J R，Delos J B，F(xiàn)lower A A，et al.Cardiovascular oscillations at the bedside：Early diagnosis of neonatal sepsis using heart rate characteristics monitoring ［J］.Physiological Measurement，2011，32 （11）：1821－1832.