相空间重构

一般的时间序列主要是在时间域中进行模型的研究,而对于混沌时间序列,无论是混沌不变量的计算,混沌模型的建立和预测都是在所谓的相空间中进行,因此相空间重构就是混沌时间序列处理中非常重要的一个步骤。

所谓混沌序列,可以看作是考察混沌系统所得到的一组随着时间而变化的观察值。假设时间序列是 \(\{ x(i): i=1,\cdot \cdot \cdot, n\}\),那么吸引子的结构特性就包含在这个时间序列之中。为了从时间序列中提取出更多有用的信息,1980年Packard等人提出了时间序列重构相空间的两种方法:导数重构法和坐标延迟重构法。而后者的本质则是通过一维的时间序列\(\{x(i)\}\)的不同延迟时间 \(\tau\) 来建\(d\)维的相空间矢量

\[y(i)=(x(i),...,x(i+(d-1)\tau)), 1\leq i\leq n-(d-1)\tau\]

1981年Takens提出嵌入定理:对于无限长,无噪声的\(d\)维混沌吸引子的一维标量时间序列\(\{x(i): 1\leq i \leq n\}\)都可以在拓扑不变的意义下找到一个\(d\)维的嵌入相空间,只要维数\(d\geq 2d^{'}+1\)。根据Takens嵌入定理,我们可以从一维混沌时间序列中重构一个与原动力系统在拓扑意义下一样的相空间,混沌时间序列的判定,分析和预测都是在这个重构的相空间中进行的,因此相空间的重构就是混沌时间序列研究的关键。

I. 嵌入维数和延迟时间

参考文献:时间序列模型之相空间重构模型

关于嵌入维数\(d\)和延迟时间\(\tau\)的取值,通常有两种观点:第一种观点认为\(d\)\(\tau\)是互不相关的,先求出延迟时间之后再根据它求出合适的嵌入维数。求出延迟时间\(\tau\)比较常用的方法主要有自相关法,平均位移法,复自相关法和互信息法等,关键的地方就是使得原时间序列经过时间延迟之后可以作为独立的坐标来使用。同时寻找嵌入维数的方法主要是几何不变量方法,虚假最临近法(False Nearest Neighbors)和它改进之后的Cao方法等。第二种观点则是认为延迟时间和嵌入维数是相关的。1996年Kugiumtzis提出的时间窗长度则是综合考虑两者的重要参数。1999年,Kim等人提出了C-C方法,该方法使用关联积分同时估计出延迟时间和时间窗。

延迟时间的确定:自相关系数法,交互信息法

嵌入维数的确定:几何不变量法,虚假最临近点法,虚假最临近点法的改进Cao方法

I.I. C-C

1
2
%% 计算嵌入维数 Num_m 和延迟时间 Num_tau
[Num_Smean,Num_Sdeltmean,Num_Scor,Num_tau,Num_m]=C_CMethod_inf(Y_series,20);

II. 构造相空间矩阵

II.I. 方法1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
% 默认tau=1

x=data1;
learnlength=800;% 已知数据集长度
trainlength=700;% 训练数据集长度
testlength=learnlength-trainlength;% 测试数据集长度
xlearn=x(1:learnlength);

%% 训练step_to_predict步预测的模型
xtrain=x(1:trainlength);
% 对于一维数据,必须列向量
step_to_predict=10;
embdim=70; % embedding dimension
%构造全输入输出矩阵
X=windowize(xtrain,1:embdim+step_to_predict);
%inputs of the training data,取embedding dimension维数
Xtrain = X(:,1:embdim);
%outputs of the training data,取最后一列
Ytrain= X(:,end);
% 1,...,70---->80
% 2,...,71---->81
% 629,...,690---->700

%% 测试
% 能够预测的个数
ntest=testlength-step_to_predict+1;
for j=1:1:ntest
% 从trainlength开始向前预测
test_matrix=xlearn(trainlength-embdim+j:trainlength+j-1);
% 691,...,700---->710
% 692,...,701---->711
% 790,...,799---->800
Xt=test_matrix';
Yt=simlssvm({Xtrain,Ytrain,type,gam,sig2,'RBF_kernel','preprocess'},{alpha,b},Xt);
Ytest(j)=Yt;
end

II.II. 方法2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
function [xn,dn] = PhaSpaRecon(s,tau,m)

% 混沌序列的相空间重构 (phase space reconstruction)
% [xn,dn] = PhaSpaRecon(s,tau,m)
% 输入参数: s 混沌序列
% tau 重构时延
% m 重构维数
% 输出参数: xn 相空间中的点序列(每一列为相空间中一个点)
% dn 一步预测的目标

len = length(s);
if (len-1-(m-1)*tau < 1)
disp('err: delay time or the embedding dimension is too large!')
xn = [];
dn = [];
else
xn = zeros(m,len-1-(m-1)*tau);
for i = 1:m
xn(i,:) = s(1+(i-1)*tau : len-1-(m-i)*tau); % 相空间重构,每一列为一个点
end
dn = s(2+(m-1)*tau : end); % 预测的目标

% 举例分析
% xn的列向量作为输入,dn的列向量作为输出
% 相空间重构,每一列为一个点
% len=100,tau=1,m=10
% xn(1,:)=s(1:90)
% xn(2,:)=s(2:91)
% xn(10,:)=s(10:99)
% dn=s(11:end)
% 1,2,...,10---->11
% 2,3,...,11---->12
% 90,91,...,99---->100

% len=100,tau=3,m=10
% xn(1,:)=s(1:72)
% xn(2,:)=s(4:75)
% xn(10,:)=s(28:99)
% dn=s(29:end)
% 1,4,...,28---->29
% 2,5...,29---->30
% 72,75,...,99---->100

end

构造训练数据集、测试数据集。

1
2
3
4
5
interval = 4;
[xn_train,dn_train,xn_test,dn_test] = TrainTestSample(interval,xn,dn);

% 相空间中点的轨迹分解成:训练样本与测试样本
% 在所有样本中,每隔 interval 个样本为训练样本,其余为测试样本
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
function [TrainFeature,TrainKnnTarget,...
TestFeature, TestKnnTarget]...
= TrainTestSample(interval,FeatureArray,KnnTarget);
% 将特征矩阵分解成:训练样本与测试样本

[row,column] = size(FeatureArray);
TrainColumn = 1:interval:column; % 训练样本的列数
Column = 1:column;
Column(TrainColumn) = []; % 删除训练样本的列数,余下的是测试样本列数
TestColumn = Column; % 测试样本的列数

%========= 训练样本 ========%

TrainFeature = FeatureArray(:,TrainColumn); % 训练样本
TrainKnnTarget = KnnTarget(:,TrainColumn); % 训练目标

%========= 测试样本 ========%

TestFeature = FeatureArray(:,TestColumn); % 测试样本
TestKnnTarget = KnnTarget(:,TestColumn); % 测试目标

II.III. 方法3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function Data=reconstitution(data,N,m,tau)
%该函数用来重构相空间
% m为嵌入空间维数
% tau为时间延迟
% data为输入时间序列
% N为时间序列长度
% Data为输出,是m*n维矩阵
M=N-(m-1)*tau;%相空间中点的个数
Data=zeros(m,M);
for j=1:M
for i=1:m %相空间重构
%Data(i,:)=data(((i-1)*tau+1):1:((i-1)*tau+M));
Data(i,j)=data((i-1)*tau+j);
end
end