简介：相空间重构是复杂系统动力学研究的基础手段，用于从单变量时间序列恢复系统完整状态，揭示混沌、周期与稳定特征。这套MATLAB资源适合学习非线性时间序列分析、需要掌握Cao法、FNN、互信息等算法的学生和科研人员，避免陷入理论公式而难以落地。压缩包共2个文件，包含1个MATLAB脚本和1个TXT说明文档，整包仅1KB，轻量便于阅读与修改。脚本覆盖数据清洗与标准化、延迟时间计算、嵌入维数确定、相空间延拓重建以及轨迹绘图等步骤，文档则进一步解释各参数含义与调参建议。基于此可快速构建相空间，并扩展计算Lyapunov指数、绘制分岔图与吸引子，用于气象、生物医学、金融等实际数据诊断。目前已有1860人学习下载，无论是入门理解还是MATLAB实验验证，都是一份简洁实用的参考。

1. 相空间重构：从一维观测序列里“还原”动力系统的道与术

做振动诊断、脑电分析或者负荷预测的人，手里最容易握着的是一根一维时间序列：传感器每秒采一个数值，采了上万点。可单看这条曲线，你只能得到“现在是多少”，得不到“系统现在处在什么状态”。相空间重构理论（Takens嵌入定理）给出的答案是：把一个观测点的历史延迟组合成向量，就能在高维空间里还原出原动力系统的吸引子结构，只要嵌入维数和延迟时间取得对，重构出的相空间与原系统拓扑等价。拿MATLAB做这件事，成本低、可视化直观，做混沌时间序列预测和故障特征提取的人都绕不开它。这篇文章从定参数到写代码，把整条链路掰开讲清楚。

2. 嵌入维数与延迟时间：Takens定理之后的两个“玄学”参数怎么定

Takens定理只断言“存在一个足够大的嵌入维数”，却没说这个维数是多少，也没提延迟时间取几拍。于是这两个参数成了整个相空间重构里最像“玄学”的部分：一组\tau和m，能把吸引子画成结构清晰的涡旋，也能画成乱麻。实际工程里，延迟时间\tau决定相邻状态向量携带多少新信息，嵌入维数m决定状态空间能展开到几维。两个参数之间有耦合，不是孤立的。下面用两种最常用的手段来定参：C-C方法（关联积分法）和虚假邻近点法（FNN）。前者擅长定延迟时间\tau，后者擅长定嵌入维数m。

2.1 C-C方法：先用关联积分把延迟时间定下来

C-C方法的核心思想是：用关联积分统计序列在不同延迟下的自相关程度，找一个“信息不再明显增加、又没过早饱和”的时刻作为延迟时间。具体做法是把长度为N的序列先去做相空间重构的“试探”——对每个候选延迟\tau，计算关联积分C(m, N, r, \tau)，再构造统计量S(\tau)和修正统计量\Delta S(\tau)。当S(\tau)第一次过零，或者\Delta S(\tau)达到第一个局部极小值时，对应的\tau就是最佳延迟。这个方法是同时估计\tau和m的，但实际用的时候，大家更常拿它先定\tau。

下面是一份不依赖任何工具箱的C-C方法实现，输入时间序列、最大延迟和候选嵌入维数范围，输出统计量和最优\tau：

function [tau_opt, S, dS] = cc_method(x, max_tau, max_m) % C-C方法：用关联积分统计量估计延迟时间tau % 输入: % x 列向量，长度N的一维观测序列 % max_tau 最大延迟搜索范围，典型取 10~20 % max_m 候选嵌入维数上限，典型取 5~10 % 输出: % tau_opt 最优延迟时间 % S 平均统计量随tau变化的向量 % dS delta S 统计量随tau变化的向量 N = length(x); % 对候选tau计算统计量 for tau = 1:max_tau % 将序列分成 tau 个互不重叠的子序列 S_sum = 0; for s = 1:tau sub_x = x(s:tau:end); % 计算该子序列的关联积分 C1 = correlation_integral(sub_x, 1, tau); for m = 2:max_m Cm = correlation_integral(sub_x, m, tau); S_sum = S_sum + (Cm - C1^m); end end S(tau) = S_sum / (tau * (max_m - 1)); end % 计算 delta S：S(tau) 在相邻 tau 间变化量的平均值 for tau = 1:max_tau-1 dS(tau) = mean(abs(S(tau+1) - S(tau))); end % 最佳tau取delta S第一个局部极小值位置 [~, locs] = findpeaks(-dS); if isempty(locs) tau_opt = 1; else tau_opt = locs(1); end end function C = correlation_integral(sub_x, m, tau) % 子函数：对给定子序列计算关联积分 % 用 Grassberger-Procaccia 关联积分的简化实现 N = length(sub_x); M = N - (m-1)*tau; if M < 2 C = 0; return; end % 构造重构向量 Y = zeros(M, m); for i = 1:m Y(:, i) = sub_x(1+(i-1)*tau : M+(i-1)*tau); end % 距离阈值 r 取数据标准差的 1/4，这是C-C方法的常见取法 r = std(sub_x) / 4; count = 0; total = M*(M-1)/2; for i = 1:M-1 dist = sqrt(sum((Y(i+1:end, :) - Y(i, :)).^2, 2)); count = count + sum(dist < r); end C = 2 * count / total; end

这段代码最关键的部分在correlation_integral子函数：它对子序列做了延迟嵌入，再统计半径r内的点对比例，得到关联积分。公式里的r取标准差四分之一是C-C方法里的经验值，不唯一，但在这个算法里足够稳定。cc_method主函数里把不同嵌入维数m的关联积分做了平均，得到S统计量，再通过\text{dS}的局部极小值找\tau。用findpeaks(-dS)是取反后找峰，正好对应“Delta S的极小值”。如果不想调findpeaks，也可以用一行[~, idx] = min(dS)直接取全局最小，但局部极小更符合算法原意。

跑这段代码时，x必须是列向量，行向量要提前转置。max_tau不要超过N/20，否则后半段子序列长度太短，关联积分失去统计意义。如果结果里\tau=1，先别急着认为序列不可重构，很可能是数据本身信噪比太低，或者子序列划分出了边缘效应。

2.2 虚假邻近点法：嵌入维数不是越大越好

嵌入维数定得过小，高维吸引子会被压扁到低维空间，原本不相邻的两个点在投影后会变成“邻居”，这些邻居叫虚假邻近点（False Nearest Neighbors）。随着m增大，虚假邻居比例会快速下降；当比例降到接近0或者不再明显变化时，那个m就是合适的嵌入维数。和C-C方法比，FNN更直接：它不需要估计任何延迟，只关心维数对空间结构的影响。

function m_opt = fnn_method(x, tau, max_m, tol) % 虚假邻近点法：估计嵌入维数m % 输入: % x 列向量时间序列 % tau 延迟时间，可由C-C方法预先得到 % max_m 最大候选嵌人维数 % tol 虚假邻居判定阈值，典型取10 % 输出: % m_opt 推荐的嵌入维数 N = length(x); fnn_ratio = zeros(max_m, 1); for m = 1:max_m % 重构m维相空间 M = N - (m-1)*tau; % 相空间点数 if M < 10 break; end Y = zeros(M, m); for i = 1:m Y(:, i) = x(1+(i-1)*tau : M+(i-1)*tau); end % 对每个点找最近邻，判断是否虚假 fnn_count = 0; for i = 1:M-1 dist = sqrt(sum((Y(i+1:end, :) - Y(i, :)).^2, 2)); [d_min, idx] = min(dist); if d_min == 0 % 距离为0的退化情况说明有重复采样，不计入 continue; end j = i + idx; % 最近邻在原始索引中的位置 if j + tau <= N && m < max_m dist_next = abs(x(i + m*tau) - x(j + m*tau)); if dist_next > tol * d_min fnn_count = fnn_count + 1; end end end fnn_ratio(m) = fnn_count / (M-1); end % 找到比例首次低于 5% 的 m，且之后不再回升 m_opt = find(fnn_ratio < 0.05, 1, 'first'); if isempty(m_opt) m_opt = max_m; end end

这个函数对每个候选m，把每个点的最近邻的下一时刻距离和当前距离做比较。如果下一时刻的距离比当前距离放大了10倍以上，说明这两个点在提升维数后可能被分开，它们就是虚假邻居。阈值tol取10是Williams的经典默认值，对多数工程序列都成立，但如果数据噪声很大，可以放宽到12或15，收严到8也行。找到fnn_ratio第一次低于5%的m，就是推荐嵌入维数。这个5%是经验门槛，取3%会偏保守，取10%更容易拿到小m。

实际使用中，FNN对\tau很敏感。如果\tau取得太大，吸引子轨道被“拉直”，最近邻本身失去意义，FNN比例可能一直降不下去；所以建议先跑C-C定\tau，再跑FNN定m。反过来，如果你已经有一个工程约定（比如振动信号常用m=5~7），FNN可以当作复核工具看这个约定是否合理。

2.3 参数联合寻优的替代路线：网格扫描与人工判读

C-C和FNN各自给出一个推荐值，但两个参数其实绑定在一起。有时\tau=5配m=4好用，\tau=3配m=6反而更好。这时候我会做网格扫描：\tau取1到15，m取2到10，对每组参数重构相空间，计算吸引子的“展开度”——比如对重构矩阵做奇异值分解（SVD），看前几个奇异值之和占全部奇异值之和的比例。比例越大，说明信息集中在少数几个主方向上，吸引子展开得越好。展开度最高的参数组合，通常就是工程上最合适的。

tau_range = 1:15; m_range = 2:10; score = zeros(length(tau_range), length(m_range)); for ti = 1:length(tau_range) for mi = 1:length(m_range) tau = tau_range(ti); m = m_range(mi); if (m-1)*tau >= length(x) % 序列长度不够 score(ti, mi) = 0; continue; end Y = reconstructPhaseSpace(x, m, tau); % 见3.4节 [~, sing, ~] = svd(Y, 'econ'); score(ti, mi) = sum(sing(1:3)) / sum(sing); end end [row, col] = find(score == max(score(:))); tau_best = tau_range(row(1)); m_best = m_range(col(1));

这个网格扫描本质上是把相空间重构当成一个“特征提取器”在评估，SVD的奇异值衰减越慢，说明吸引子结构越复杂，越值得用更高维去展开。如果score的值很接近（差不到5%），优先取更小的m，因为维数低意味着后续预测模型参数更少，更不容易过拟合。这和FNN的结论通常一致，但当序列长度有限时，网格扫描比FNN更稳健——FNN在高维下需要足够的近邻点，序列短时统计量会失真。

2.4 定参之后的可视化：吸引子重构质量怎么看

参数定完，第一件事不是急着建模，而是画图。把三维相空间轨迹画出来，好吸引子和坏吸引子一眼能分辨：好的重构会呈现光滑的流形结构——比如Lorenz系统的双螺旋、Rossler系统的单环，轨道不会交叉成一团；差的重构要么是一条贴在对角线附近的细带（\tau太小，信息冗余），要么是铺满整个空间的噪声球（\tau过大或m过高）。我习惯先用plot3看三维视图，再用scatter看二维投影，两张图都过目了才继续往下走。这里不截图，但你完全可以拿一段实际振动数据试，效果非常直观。

3. 用MATLAB把重构流程接起来：完整代码与每行注释

参数理论说了一堆，落地才是硬道理。这一章给一条完整的MATLAB流程，从读数据到最后画相图，每个环节都能直接复制到你的脚本里。

3.1 读入时间序列与数据预处理：去趋势、去噪、归一化

相空间重构最容易被轻视的环节是预处理。原始传感器数据常常带线性漂移、直流偏置和异常尖峰。漂移会让C-C方法的关联积分计算失真——序列缓慢上升时，距离阈值r被抬高，虚假邻近点判定错乱。所以第一步要detrend去掉线性趋势，再减均值、除标准差做z-score归一化。归一化不是必须的，但它能让\tau和m的搜索结果对不同量纲的序列有一致性，尤其是C-C方法里r取标准差比例，量纲变化会影响关联积分的绝对数值。

data = readmatrix('sensor_data.csv'); % 读入数值矩阵 t = data(:, 1); % 第一列通常是时间戳 x = data(:, 2); % 第二列是观测值，取列向量 x = detrend(x, 'linear'); % 去线性趋势，消除漂移 x = x - mean(x); % 去掉均值，让序列围绕0波动 x = x / std(x); % 归一化到单位方差

预处理不是越猛越好。detrend只做线性趋势去除，别用高阶多项式去拟合趋势，那会把低频动力学信息也当趋势去掉。去噪用移动平均要谨慎，窗口大了会把吸引子的精细结构抹平；如果数据是轴承振动这类高频信号，我一般最多做一次5点移动平均，或者干脆不去噪，让后续的相空间重构自己去处理噪声。归一化这步不改变相空间重构的数学本质，但它让cc_method里的r取std(x)/4变得有明确含义，也让不同段数据的参数可比。

3.2 C-C方法计算延迟时间：一个可复现的m函数

把第2.1节的cc_method存成独立m文件保存到当前目录，然后在主脚本里调用：

max_tau = 15; % 最大延迟：别超过序列长度的1/20 max_m = 6; % 候选嵌入维数上限 [tau_opt, S, dS] = cc_method(x, max_tau, max_m); % 画统计量曲线，人工复核最优tau figure; subplot(2,1,1); plot(1:max_tau, S, 'o-'); grid on; xlabel('\tau'); ylabel('S(\tau)'); title('S统计量'); subplot(2,1,2); plot(1:max_tau-1, dS, 'o-'); grid on; xlabel('\tau'); ylabel('\Delta S(\tau)'); title('Delta S统计量');

这段代码把C-C方法的输出直接画出来。S(\tau)曲线出现明显的第一次过零位置，或者\Delta S的第一个尖峰位置，就是\ta_{\mathrm{opt}}。自动取局部极小时，偶尔会取到\tau很小但并非全局最优的点，所以画图人工复核是习惯动作。如果发现自动结果和曲线肉眼看到的极值不一致，以曲线为准，手动在下一行代码里覆盖tau_opt就行。

3.3 虚假邻近点法确定嵌入维数m

有了延迟时间，接下来定嵌入维数。把2.2节的fnn_method也存成m文件，调用如下：

max_m = 10; tol = 10; % 距离放大倍率阈值 m_opt = fnn_method(x, tau_opt, max_m, tol); fprintf('C-C得到 tau_opt = %d\n', tau_opt); fprintf('FNN得到 m_opt = %d\n', m_opt);

FNN的结果如果恰好是max_m边界值（比如你设了10，它返回10），说明可能在当前\tau下虚假邻居比例还没收敛到5%以下，这时候把max_m提高到15再跑一次。FNN依赖最近邻搜索，序列长度越短，最近邻的可靠性越差；要求M至少大于50，也就是N - (m-1)*tau要大于50，否则低维段的结果都不可信。

3.4 相空间重构主函数：从状态向量到嵌入矩阵

参数定了，下一步就是把一维序列扩展成多维状态矩阵。常见做法是自己写一个重构函数，不依赖特定工具箱，逻辑清晰也好改。下面这个函数我用了很多年，核心就是按列取延迟：

function Y = reconstructPhaseSpace(x, m, tau) % 相空间重构：把一维序列扩展为m维延迟嵌入矩阵 % 输入: % x N×1列向量，预处理后的观测序列 % m 嵌入维数 % tau 延迟时间 % 输出: % Y (N-(m-1)*tau) × m 矩阵，每一行是一个状态向量 N = length(x); nPoints = N - (m-1)*tau; if nPoints < 1 error('序列长度不足以重构：需要至少 (m-1)*tau+1 个点'); end Y = zeros(nPoints, m); for i = 1:m % 第 i 列 = 原序列从第 (i-1)*tau+1 个点开始取 nPoints 个点 Y(:, i) = x((i-1)*tau+1 : (i-1)*tau+nPoints); end end

逻辑很简单，但有个关键点容易搞错：nPoints是重构后的相空间点数，每一列按延迟递增偏移(i-1)*tau。比如原始序列有1000个点，m=3，tau=5，那么输出矩阵是982行×3列，丢掉的是开头和结尾各一段无法完整对齐的点。丢掉头尾是很正常的，预测时最后(m-1)*tau个点的“未来值”不存在，只能用来补全输入而不能作为训练目标。

3.5 重构后的吸引子三维轨迹图与相图矩阵

重构完成后立刻画图。三维吸引子用plot3，前面两条坐标是重构矩阵的第一列和第二列，竖直方向取第三列。如果你m大于3，可以用subplot矩阵画不同列组合的二维相图：

Y = reconstructPhaseSpace(x, m_opt, tau_opt); figure; plot3(Y(:,1), Y(:,2), Y(:,3), 'b-', 'LineWidth', 0.3); xlabel('x(t)'); ylabel('x(t+\tau)'); zlabel('x(t+2\tau)'); title(sprintf('3D相空间吸引子 (m=%d, tau=%d)', m_opt, tau_opt)); grid on; axis tight; % 如果m>3，画前四列的二维组合 if m_opt >= 4 figure; for p = 1:4 subplot(2,2,p); plot(Y(:,1), Y(:,p+1), '.', 'MarkerSize', 1); xlabel('x(t)'); ylabel(sprintf('x(t+%d\\tau)', p)); grid on; end end

如果吸引子是成片的密集点而不是连续线条，可以把线型改成点，或者把MarkerSize调小。连续光滑的轨道说明数据采样率足够高，动力学信息完整；如果画出来像散沙，别急着怀疑代码，先回头检查\tau是不是被C-C方法算成了一个过大的值。把\tau调小一至两拍再画，通常会看到结构。

4. 相空间重构避坑指南：五条能让你重构结果直接翻车的常见问题

做相空间重构最大的错觉是“参数算出来就能直接造模型”。实际上我从第一次调C-C方法到现在，踩过的坑没有一个是可以跳过的。下面五条按“现象→原因→解决”写清楚。

4.1 延迟时间过大导致吸引子“拉平”，过小导致信息冗余

现象：\tau取大了，三维相图变成一条几乎贴近主对角线的细带，或是一团摊开的云；\tau取小了，相图接近一条直线，轨道自相交严重。 原因：\tau太大，相邻状态向量之间的相关性完全丢失，每一步延迟相当于从噪声中取样，吸引子被“拉散”；\tau太小，状态向量的各维几乎相同，等于反复使用同一个信息点，空间没有展开。 解决：以C-C方法的结果为中心，左右各扫描两三个值，对比相图结构。我一般取\tau作为初值，再测试\tau±1、\tau±2，选择能把吸引子展开成最清晰轨道的那一个。别只用自动结果，C-C的极值点和可视化的最优值经常相差一拍。

4.2 嵌入维数过高导致噪声被放大，计算量陡增

现象：FNN推荐m=3，你为了“捕捉更多信息”改成m=6，结果预测模型训练时间翻了倍，准确率却掉了5个点。 原因：嵌入维数高会把噪声也嵌入进去。每增加一维，状态空间体积指数级增长，需要的数据量也指数级增长；数据量不足时，高维空间里大部分区域是空的，近邻关系丧失物理意义。 解决：FNN给出m后，顺带算一下数据量是否够：要求相空间点数nPoints至少是m的20倍。不满足就降维，或者去采更多数据。训练预测模型时，从m_opt开始只做减法，不做加法。

4.3 数据未预处理：趋势项和异常值让C-C方法失效

现象：同样一段信号，原始数据跑C-C方法得到\tau=1，去掉趋势和归一化后得到\tau=7，差距巨大。 原因：C-C方法依赖关联积分统计点对距离，线性趋势会让距离分布严重偏斜，r取固定比例的分位数时，大部分点对距离落在r之外，关联积分失去分辨力。异常尖峰也会产生同样的效应。 解决：重构前必须做detrend和归一化，但别做过度滤波。尤其不要用滑动平均窗口超过20点的滤波器，那会把吸引子的高频细节全部抹掉。

4.4 重构后的数据尺度差异导致后续模型训练翻车

现象：相空间重构矩阵各列的实际数值范围完全一样，但和别的特征拼接（比如均值、峰峰值）之后，尺度过大的列主导了模型损失，导致训练不收敛。 原因：重构矩阵各列来自同一个序列，天然同尺度，但如果你把重构结果作为特征输入LSTM或神经网络，往往还要拼接其他统计特征；后者量纲差异会打乱梯度更新。 解决：重构矩阵先保持原始尺度用于可视化，拼特征时单独做一次zscore标准化。别把重构矩阵自己再归一化，不然吸引子结构的相对距离会被破坏，影响KNN这类依赖距离的模型。

4.5 官方函数与自写函数的边界：什么时候别用自带实现

现象：MATLAB自带的phaseSpaceReconstruction在某些版本里返回的结果列顺序或行数与你手写的延迟矩阵不同，导致下游代码报错。 原因：官方函数更侧重于可视化场景，返回的是降维约简后的相空间轨迹，不是单纯的延迟嵌入矩阵，行数和延迟对齐方式有自己的约定。 解决：如果只是想快速看一眼吸引子，用自带的；如果要拿相空间矩阵去做训练集构造，一定用自己的reconstructPhaseSpace，逻辑透明、可复现。两者结果不等价不是bug，是定位不同。我还见过有同学把官方函数输出直接当训练输入，结果序列长度对不上，报错后才换回来。工程上我主张：建模阶段永远自写重构函数，画图阶段可以随意。

5. 从重构到落地：拿重构后的相空间去做预测与故障诊断的验证技巧

重构得对不对，最终要看下游任务是否受益。我最常用的验证手段是做一个“重构+KNN单步预测”的基线实验。流程很简单：先把序列重构为Y矩阵，然后把Y的前n-1行作为输入，把原始序列中对应的“下一时刻观测值”作为标签，训练一个KNN回归模型，看测试集上的RMSE。如果重构参数选得好，这个基线的RMSE会明显低于直接用原始一维序列做的滞后预测。

% 用重构相空间的前 nPoints-1 行做训练，预测下一步观测 Y = reconstructPhaseSpace(x, m_opt, tau_opt); nPoints = size(Y, 1); % 目标值：每个相点对应的下一个观测值 target_start = (m_opt-1)*tau_opt + 2; % 原序列中与Y第一行对齐的下一点 target = x(target_start : target_start + nPoints - 1); X_train = Y(1:end-1, :); y_train = target(1:end-1); X_test = Y(end, :); % 用最新状态预测后一步 y_true = target(end); mdl = fitcknn(X_train, y_train, 'NumNeighbors', 5, 'Standardize', 1); y_pred = predict(mdl, X_test); fprintf('KNN单步预测RMSE = %.4f\n', sqrt(mean((y_pred - y_true).^2)));

这是一个基线验证，不是最终预测方案。KNN对距离度量敏感，所以Standardize要开；邻居数小数据默认5即可，数据量大可以调到10。如果你在做LSTM或Transformer，重构矩阵直接作为时间窗口输入也会有提升——很多从业者看到这里才明白，LSTM的“隐状态”本质上也是在学一个非线性嵌入，只是参数不可解释，相空间重构给了它一个可解释的初始化视角。

故障诊断场景里，另一个验证技巧是看递归图（Recurrence Plot）。把重构后的状态向量两两计算距离，画成灰度图，健康状态的递归图通常有规律的纹理，故障状态会出现明显的色块断裂或带状结构。这是相空间重构在故障诊断里最讨巧的一种应用：不需要训练任何模型，一眼就能看出状态是否发生突变。做法也简单，pdist2(Y,Y)之后用imagesc显示阈值化结果。

我自己最早做轴承故障诊断时，偷懒直接用\tau=1、m=2，画出来的相空间挤成一个圆饼，怎么看都像高斯噪声。后来老老实实跑C-C和FNN，再把递归图打开，才真正看到不同故障类型在相空间轨迹上的分岔差异。每次换新数据，我都会把这张递归图打印出来贴在工作笔记上，确认重构质量过关了再继续建模。参数自动计算只是起点，花十分钟画图做人工复核，是那一整套流程里最值回票价的一步——希望帮到你。

本文还有配套的精品资源，点击获取