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

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

0 引 言

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

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

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

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

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

1.2 樣本熵

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

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

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

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

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

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

對于m＋1點矢量，同樣有

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

步驟4 此序列的樣本熵為

但是在實際計算中N 不可能為∞，當N 取有限值時，估計

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

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

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

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

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

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

Lorenz方程公式如下所示

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

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

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

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

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

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

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

3 實驗結(jié)果及分析

3.1 Lorenz

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

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

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

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

3.2 Logistic離散映射

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

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

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

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

3.3 一階自回歸AR 模型

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

同前兩個實驗一樣，本次也計算了AR 模型產(chǎn)生的8種長度數(shù)據(jù)的特征量，圖3 （a）、圖3 （b）僅給出數(shù)據(jù)長度N 分別為300和800時第1次產(chǎn)生的原始數(shù)據(jù)和對應的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時的t檢驗臨界值），表明原始數(shù)據(jù)與替代數(shù)據(jù)近似熵與樣本熵的值的差異沒有達到統(tǒng)計檢驗顯著性，即原始數(shù)據(jù)和替代數(shù)據(jù)不存在顯著差異，可以判斷該時間序列中不包含非線性成分，是線性的 （和AR 時間序列已知特性相符），并且該結(jié)果不受數(shù)據(jù)長度的影響，另外樣本熵計算用時比近似熵要短。

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

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

4 結(jié)束語

論文在三組仿真數(shù)據(jù)上，驗證了樣本熵是否可以作為特征量來驗證時間序列中有無非線性成分。實驗結(jié)果表明，對于不同長度的時間序列，樣本熵均可以作為替代數(shù)據(jù)法的特征量檢驗出時間序列中是否含有非線性成分；并且采用樣本熵的替代數(shù)據(jù)法比采用近似熵的替代數(shù)據(jù)法計算效率得到明顯提高。所以，本研究得到結(jié)論：基于樣本熵的替代數(shù)據(jù)法是一種準確、高效的非線性檢測方法。當然本研究還需要進一步在其它數(shù)據(jù)上進行驗證，來支持論文的結(jié)論。但是本文的結(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］.河西學院學報，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］.長春師范學院學報（自然科學報），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）.［于青.關聯(lián)維數(shù)計算的分析研究［J］.天津理工學院學報，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.