非线性动力学数据分析

时间序列分析读书报告与数据分析

刘愉

200921210001

时间序列分析是利用观测数据建模，揭示系统规律，预测系统演化的方法。根据系统是否线性，时间序列分析的方法可分为线性时间序列分析和非线性时间序列分析。

一、时间序列分析涉及的基本概念

对于一个动力系统，我们可以用方程表示其对应的模型，如有限差分方程、微分方

1、测量

程等。如果用Xt或X(t)表示所关心系统变量的列向量，则系统的变化规律可表示成

Xt1f(Xt)或

dXdtF(X)

其中X可以是单变量，也可以是向量，F是函数向量。通过这类方程，我们可以研究系统的演化，如固定点、周期、混沌等。

在实际研究中，很多时候并不确定研究对象数据何种模型，我们得到的是某类模型（用Xt或X(t)表示）的若干观测值（用Dt或D(t)表示），构成观测的某个时间序列，我们要做的是根据一系列观测的数据，探索系统的演化规律，预测未来时间的数据或系统状态。

2、噪声

测量值和系统真实值之间不可避免的存在一些误差，称为测量误差。其来源主要有三个方面：系统偏差（测量过程中的偏差，如指标定义是否准确反映了关心的变量）、测量误差（测量过程中数据的随机波动）和动态噪音（外界的干扰等）。 高斯白噪声是一类非常常见且经典的噪声。所谓白噪声是指任意时刻的噪声水平完全独立于其他时刻噪声。高斯白噪声即分布服从高斯分布的白噪声。这类噪声实际体现了观测数据在理论值（或真实值）周围的随机游走，它可以被如下概率分布刻画：

p(x)dx1222exp(xM)22dx

（1）

其中M和均为常数，分别代表均值和标准差。

3、均值和标准差

最简单常用的描述时间序列的方法是用均值和标准差表示序列的整体水平和波动情况。 （1）均值

如果M是系统真实的平均水平，我们用观测的时间序列估计M的真实水平方法是：认为N个采样值的水平是系统水平的真实反映，那么最能代表这些观测值（离所有观测值最近）的Mest即可作为M的估计。于是定义Dt与Mest的偏离为(DtMest),所以，使下面E最小的M的估计值即为所求：

N22E(Dtt1Mest)

（2）

1/11 经过求道计算，得到

M1NNestDtt

1（3）

即样本的均值即为系统真是均值的估计值。

（2）标准差

标准差代表了系统在均值两侧的波动情况。对时间样本有：

VtDtMest

（4）

为了分析所有时间上平均的波动情况，我们也可以尝试对波动取平均，即：

1NNt1(DtMest1)NNt1DtMest0

（5） 我们发现，这样平均的结果是正负波动抵消了，波动的平均恒为零，为了避免这种情况，改用波动的平方的平均水平代替，即

21NNVtt121NN(Mestt1Dt)

（6）

2即为标准差。 （3）均值的标准误差

我们用Mest估计M，存在一定偏差或不确定性，即：

MestMuncertainty

（7）

实际上，这种不确定性来自每次测量偏差的平均，通常每次测量偏差是服从高斯分布的，所以平均的不确定性计算得：

N

（8）

我们称之为均值的标准误差。

二、线性时间序列分析方法及模型举例

对于线性时间序列，主要的分析方法有：均值和标准差、线性相关分析和功率谱分析。

1、均值和标准差分析前面已经讲过；

例：模型一（模型本身是确定的（无外界干扰等随机波动），观测序列是真实值加上高斯白噪声；）

有限差分方程系统：xt1Axt，其平稳状态为xtA/(1)M；观测时间序列DtxtWt，其中，Wt 独立的服从均值为0，标准差为的高斯分布。从系统的差分方程我们可以看到，系统本身不受外界干扰，是确定性模型。所以观测得到时间序列的波动完全来自于测量过程。

对于上述模型，可以通过均值、方差的估计即可估计模型、作出预测。

2、线性相关分析

2/11 这种分析方法用于研究时间上相关的序列，即后一时刻的值完全或部分由前一时刻的或前几个时刻的值决定。在模型一中，我们假设Wt之间是独立的；当这种假设不成立时，取另一种极端，即后一时刻完全取决于前一时刻的值：

Vt1f(Vt)

（9）

我们以简单的线性函数为例：

Vt1Vt

（10）

如果结合完全独立的情形与式（10），则有以下情况：

Vt1VtW

（11）

ρ在-1到1之间取值，ρ越接近0，数据间越不相关；ρ接近1，表示线性正相关；ρ接近-1，表示线性负相关

通过时间序列的一系列观测值Dt减去均值得到Vt，我们可以通过以下公式计算相关系数，

estt1N1Vt1VtVtVtt1N

1（12）

例：模型二（模型本身有不确定因素（外界干扰），观测序列是真实值加上高斯白噪声）

受外界因素影响的有限差分方程：xt1Axtvt，引入的vt是外界干扰造成的系统本身的波动，测量过程仍然像Model One一样，DtxtWt，这是如果做Vt1对Vt的变化图（见课本figure 6.7），发现二者之间有强烈的线性关系。对于这类模型，我们即可用线性相关分析来建模、预测。

如果将线性相关加以推广，可以得到自相关函数，它反映的是Vt与Vtk之间的关系：

R(k)t1NkVtkVtVtVtt1Nk

（13）

3、功率谱分析 （1）傅里叶变换

对线性系统，一个信号可以分解成为不同频率的正弦波。

（a）频率为ω的正弦输入，它的输出也是同频率的正弦信号，但是幅度和相位可能发生改变。输出正弦波的振幅与输入正弦波的振幅满足：

Aoutput()G()Ainput()

（14）

输出相位相对输入相位在每个频率上有固定的偏移，即：

()output()input()3/11

（15） G()称为系统的增益，它在不同频率上通常不一样。()称为相移，在不同的频率成分通常相移也不同。

（b）线性叠加的输入的输出结果等于各个输入分别输入时的输出的叠加。把一个信号分解成不同频率正弦信号的方法即傅里叶变换。

特殊的，输入为白噪声时，Ainput()是一个与噪声标准差成正比的常数，与频率无关，即白噪声可以认为是所有不同频率成分信号之和，所以称之为“白”。 （c）传输函数

如果已知输入和输出，可以得到：

G()Aoutput()Ainput()

()outp()inpu()utt

（16）

（16）中两个函数成为出传输函数，可以用于描述系统特性。

（d）功率谱

如果我们不能准确得到输入信号，但是我们知道或假设它是白噪声。则Ainput()就是常数，进而有：

G()constAoutput()

（17）

G()的平方称为功率谱。功率谱包含了与自相关函数完全一样的信息。事实上，功率谱就是自相关函数的傅里叶变换。尽管它们蕴含的信息是一样的，但不同形式使它们在分析数据时又具有各自的优势。所以有时使用功率谱来分析数据比用自相关函数更有优势。

三、非线性时间序列分析方法及模型

前面列举了一系列线性时间序列的分析方法，但是对于非线性系统，存在一种特殊的状态，即混沌状态，对于混沌状态的时间序列，我们无法用线性的分析方法区分。

例：第一章的有限差分方程：

xt1xt(1xt)

（18）

Dtxt

（19）

观测值即使不引入噪声，其时间序列也在不断波动，当=4时，系统进入混沌状态，用线性自相关函数分析，如图6.14，发现我们无法区分这个非线性模型与模型一。

我们需要探索一些分析非线性时间序列的方法。对于非线性时间序列分析，主要包括两部：重构系统动态模型和系统特征的刻画。

1、系统动态模型的重构：

（1） 对于有限差分方程——构建return map

Return map 是观测值Dt1关于Dt的图像（回归曲线），反映的是xt1与xt的关系。 （2） 对于微分方程——重构相平面

一维高阶微分方程可以转化为多维一阶微分方程组，以二阶微分方程为例，

4/11

dxdt22bx

（20）

转化为两个一阶微分方程组：

dxydtdybxdt

（21）

要做变量x与y的相平面，首先要做如下离散化和近似： 观测值D0,D1,，将x关于t的导数近似为：

dx(t)dtx(th)x(t)hdDtdtDthDthlimh0

（22）

其中h只能取整数，最小取1，事实上，h去较大值也可以得到合适得结果。重建相平面实际是做DthDthdxdt关于Dt的图，有时候，只可以只做Dth关于Dt的图。

dydt对于更一般的微分方程：

f(x,y),g(x,y)

（23）

虽然情况更复杂，但也可以通过这种方法重建相平面，图6.17-6.19可以说明这一过程的合理性。

（3） 嵌入时间序列

对于更高维的时间序列（p维），需要用嵌入时间序列的方法构建相平面（相空间），p维的的嵌入时间序列构成如下：

Dt(Dt,Dth,Dt2h,,Dt(p1)h)

（24）

其中p是嵌入维数，h是嵌入延迟。

经过上述三种方法，可以基本得到模型的基本特征。

2、系统特征的描述：

在模型重构后，可以通过拟合等方式对系统特征做进一步刻画。

四、混沌时间序列的刻画

混沌定义：bounded, deterministic dynamics that are aperiodic and display sensitive dependence on initial conditions.根据定义中体现的混沌系统的特征，用时间序列分析的方法研究。

1、有界

有界的定义是当时间趋于无穷时，系统永远保持在有限空间内运动。这个定义直接用于时间序列分析并不是很有效，因为测量的时间序列是时间上是有限的，变化范围也是一定的。

在时间序列分析中可以用另一种方法研究系统是否有界——稳态，即时间序列在演化过程中是否体现了相同的行为特性。相似的行为可以用均值和方差衡量。一种常用的衡量方法是将时间序列等分（三分、四分或十分等等），计算每段的均值和方差是否相近或统

5/11 计意义上可以认为相同。

如果一个时间序列是分平稳的，我们可以通过对时间序列做一些变换使之平稳，如一阶差分或后一时刻与前一时刻相除等。

2、非周期

混沌的系统是非周期的，由于噪声因素，即使周期的序列也可能出现非周期，那么如何判定时间序列是否存在周期呢？对于一维时间序列或p维的嵌入时间序列 Dt，我们定义时间i和j的测量值之间的距离定义为：

i,j|DiDj|

（25）

严格的周期T定义是当|ij|nT,n0,1,2,时，i,j0，对于有噪声的时间序列，我们定义一个距离r，当i,jr时，我们就在坐标(i, j)处打点，我们将这样做出的图成为recurrence plots，它可以看出重建的轨线如何重复自身的演化。（图6.26和图6.27是r取不同值时的图，都可以看出系统的周期，图6.28和6.29是混沌的情形）。对于混沌的时间序列，图的形状可能和r的选取有关，于是定义出这样一个correlation integral：

C(r)numberoftimes|DiDj|rN(N1)

（26）

这是一个对于混沌系统很重要的指标，它的重要意义不在于某个r处C(r)的取值，而在于C(r)如何随r变化。

（1）对周期序列，r微小的变化不会引起C(r)明显的变化；

（2）对混沌序列，r微小的变化会使C(r)明显增大，即打点明显增多； （3）对于白噪声，r微小变化时C(r)增大更快。

事实上，C(r)与分形维数密切相关，取一点做参考点，随着r增加，距离参考点r范围内的点与rv成正比，其中v是系统，所以有：

C(r)Arv

（27）

A是比例常数，两边取log得到：

logC(r)vlogrlogA

（28）

所以，只要对log C(r)与log r拟合，即可推算出维数。用相关维数可以分析混沌时间序列的吸引子，当嵌入时间维数p≥2v+1时，可以重构出系统吸引子，有时候p≥v也足够了。

3、确定性

如果已知t时刻的值，在预测下一时刻值过程中没有随机因素，那么系统就是确定的，如模型一；如果混入了随机因素（外界干扰），则系统是不确定的。但是，在观测时间序列中，噪声是不可避免的，如果预测是完美的，就成系统完全确定，如果预测是好的但不完美，就说系统有一个确定成分。

假设观测数据最后时刻是T，我们可以用一下方法预测T+1时刻的值： （1）产生嵌入时间序列Dt； （2）找到时间T的嵌入点列：

DT(DT,DTh,,DT(p1)h)

6/11 找到其他的嵌入时间序列中与DT最接近的点Da；

（3）基于系统的确定性，Da+1可以看做是由Da预测出来的，所以将Da+1作为是T+1时刻的预测值，记PT+1。

另一种预测方法是用与DT最接近的K个点Dai的下一时刻Dai1的平均值作为T+1时刻的预测值，即：

T11KKDa

1 （29）

ii1既然是预测，一定有预测误差，衡量预测误差，通常是将观测数据分为两半，用前一半数据预测后一半数据，后一半数据测量值与预测值比较，衡量预测误差，即：

1TT(Di1TkPTk)

2 （30）

ε越小，说明预测越好，至于ε小到什么程度算预测足够好，可以与最差的预测（如所有观测值序列的平均值，丧失了所有时间信息，只保留了系统的平均水平）的误差做比较，

lazy1TT(DTkk1Plazy)

（31）

当PlazyMest时，lazy其实就是时间序列的方差σ，，所以用

22即可衡量预测误差的大小，比值越小越好。

4、对初值的敏感性

混沌系统的另一重要特征就是对初值敏感，衡量系统对初值的敏感性可以用Lyapunov指数，其计算步骤可以如下概括：

（1）在给定初值后，经迭代产生序列x0,x1,,xn1 （2）计算每点处的斜率 （3）计算李氏指数：

1n1t0dfnxt dx（4） 当给定两个初始值x0，y0时，n步迭代后的值xn和yn的差距约为：

n1xnynt0dfdxxt|x0y0|

五、混沌和非线性的检测

混沌是一种复杂的现象，判定一个时间序列是否来自对非线性系统或混沌系统的观测，更严格的方法是假设检验。

1、零假设：数据来自线性系统

7/11 xt1a0xta1xt1a2xt2ap1xt(p1)vt

2、构造检验统计量：在零假设条件下用观测时间序列计算统计量，常用的三种统计量有：非线性系统的预测方差、李氏指数、相关维数等，当然还有一些其他的统计量也可以用来假设检验。

3、假设检验：在零假设下观测时间序列得到的检测统计量如果落在拒绝域内，则认为系统为非线性的或混沌的，否则接受原假设，认为是线性系统。

六、实例数据分析

P4.txt、TP8.txt是2个时间序列信号的数据文件，该数据的采样率是500Hz。试实现：

1、在时间轴上显示原始数据波形；

2、求每个信号的功率谱，在频率轴上显示结果，并对结果进行简单地讨论；

3、求每个信号的自相关函数，在时间轴显示结果，并对结果进行简单地讨论；

4、求2个信号的互相关函数，在时间轴显示结果，并对结果进行简单地讨论。 分析结果：

1、时间轴上数据波形：

从时间序列上看，两组数据基本维持在一个平衡水平，但是都存在尖峰，从时间序列看不出更丰富的信息，需要用其他方法进一步分析。

2、功率谱：

功率谱的计算有两种方式：

（1）计算时间序列的自相关函数，再对自相关函数做傅里叶变换的幅度谱； （2）时间序列傅里叶变换的幅度谱的平方除以点数N。 这里采用的第一种方法，结果如下图：

8/11

从两个时间序列的功率谱看，能量主要集中在了低频部分，高频部分能量分布极少，为了更清晰的看能量在低频的分布，我们截取0-20Hz部分的频率谱，如下图：

从图上可以看出，两个时间序列的低频成分中能量最大的频率大概都在1.7Hz左右，TP8的能量分布更集中，P4还有较多能量分布在0.5Hz左右。

3、自相关函数：

自相关函数反映的是t时刻与前t-k时刻记录值的关系，如果信号本身是周期的，其自相关函数保持与时间序列相同的周期，从下面两个序列的自相关函数图中，数据没有周期现象，而且相关函数随k增大很快降到0.2一下，并在-0.2和0.2之间震荡，TP8震荡的频率更高。

9/11

截取前1000个点，进一步观察：

从上图可以更好的体现出相关系数的变化，而且可以看出两个时间序列变化的一致性，TP8比P4相关程度衰减得更快。为了更好的研究两组数据的相关性，我们下面将做两组数据的相关函数。

4、两组数据的互相关函数：

10/11

截取两侧各500个点观察：

上图反映了两个时间序列数据的相关性，在k=0时，互相关函数最大，所以两个时间序列是同步的。

11/11

动力学分析方法

弹道动力学分析

数据分析

无碳小车动力学分析

大跨度桥梁非线性分析的论文

非线性理论学习心得

超声非线性（推荐）

武大非线性光学

非线性编辑心得

非线性编辑心得

《非线性动力学数据分析.doc》

将本文的Word文档下载到电脑，方便编辑。

推荐度：

点击下载文档