简介一份面向时间序列分析与非线性动力学研究的MATLAB源码包围绕相空间重构PSR实现三维重构流程适合信号处理、机器学习及生物医学等领域的开发者参考。包内完整实现了相空间重构的核心算法涵盖延时嵌入、塔肯斯定理、参数优选、距离计算与流形重构等关键环节并集成非线性特征分析如吸引子和混沌行为识别随附的洛伦兹系统示例便于对照验证。压缩包共十一个文件其中四个脚本文件承担核心计算与绘图三张图片直观展示三维重构结果另有说明文档、许可文件及阅读说明整体仅两百零五千字节轻量易用。已有183人学习下载适合需要快速上手MATLAB相空间重构的中高级用户可直接运行或改编用于自身数据研究。1. 相空间重构PSR三维轨迹图到底在干什么给一维序列“凭空”造出一个高维几何手里只有一路信号的时候比如一台旋转机械的振动加速度、一段脑电、一组汇率收盘价直接在时域里看波形能看出噪声和幅值但看不出系统当前到底在哪个状态。相空间重构Phase Space Reconstruction简称 PSR就是解决这个问题的把一维时间序列按延迟坐标重新展开成高维状态点让隐藏在波形背后的动力学结构自己“长”出来。当下最常用、也最容易上手的是做三维重构把轨迹画成三维曲线肉眼就能判断它是周期、混沌还是纯随机。标题里“所有代码_psr_三维重构_相空间_相空间重构_straightxx8_源码”这类包核心无非就是四件事定延迟时间、定嵌入维数、重构相空间、画图或导出特征。straightxx8 更像项目代号源码本身的工程结构与常规实现没有本质区别。这篇文章不把任何源码当黑匣子而是把这类项目最通用的实现路径拆开先讲透两个核心参数为什么难定再给一组能直接跑起来的 Python 代码最后告诉你最容易在哪里翻车以及怎么用已知系统验证你写的重构代码没白写。2. 重构前必须定死的两个参数延迟时间和嵌入维度的原理与选型2.1 延迟坐标嵌入为什么成立Takens 嵌入定理与三维重构的几何意义相空间重构的理论根基是 Takens 嵌入定理。它说的是对一个确定性动力系统如果你只观测到其中一个分量 x(t)那么用延迟坐标构造的向量 [x(t), x(tτ), ..., x(t(m-1)τ)]在 m 足够大的时候构成的几何结构与原系统的相空间拓扑等价。换句话说你虽然只装了一个传感器但可以通过时间延迟把系统的其他“隐藏变量”间接补回来。这里要提醒一个最容易误解的点三维重构并不等价于“系统的真实维数只有 3”。m 的选择由系统的吸引子维数决定m ≥ 2d1 是 Takens 定理给出的充分条件实际工程里常常 m 取 3 到 12 就够了。所谓“三维重构”常见做法是取 m3画 [x(t), x(tτ), x(t2τ)] 的三维轨迹如果 m 算出来是 5、7那就只能取前三个坐标或者做 PCA 投影到三维看。很多新手直接把 m 设成 3 硬画这是后面第 4 章要讲的典型坑。这套几何重构的价值在于周期信号重构后是闭合环拟周期信号是环面上的曲线混沌信号是分层折叠的“蝶形”或“涡旋”而随机噪声会填满整个相空间。所以看三维图不仅是为了“好看”更是做非线性时间序列分析的第一步定性判断。2.2 延迟时间 τ 的两种常用估计自相关法与互信息法附算法流程延迟时间 τ 决定了相邻两个坐标之间的信息重叠程度。τ 太小相邻坐标几乎就是同一个点轨迹会挤在主对角线附近看不出结构τ 太大相邻坐标变得完全独立轨迹摊成一团系统原本的几何关系被撕裂。工程上最常见的两种估计方法自相关法是最古老的方案计算原序列与延迟 τ 后序列的线性相关系数 R(τ)取 R(τ) 第一次降到初始值的 1/e 处或者第一次过零点对应的 τ。优点是计算快、稳定缺点是它只感知线性相关性对强非线性系统经常把 τ 选得过大或过小。互信息法是更被推荐的做法。它统计两个变量联合分布与边缘分布的 KL 散度衡量一个变量携带另一个变量的信息量。I(τ) 的第一个局部极小值代表原序列与延迟 τ 后的序列“最不相似”这个时候两个坐标之间既有信息延续又足够独立最适合作重构坐标。算法流程很固定把原序列 x_t 和延迟序列 x_{tτ} 分别分箱统计二维联合直方图 P(i,j)再算 I Σ P(i,j) · log(P(i,j) / (P(i) · P(j)))。实际处理时我一般会同时把互信息曲线画出来而不是只看程序自动挑的那个极小值。τ 的取值范围通常试探到信号特征周期的三分之一或二分之一长度采样率越高τ 的可选范围越宽。2.3 嵌入维数 m 的两种常用估计伪近邻法与 Cao 方法定完 τ 再定 m。物理直觉是这样的m 太小轨迹在高维空间本来相隔很远的点会被“压扁”到低维空间里变成近邻造成虚假的折叠m 太大噪声会被当成有效维度轨迹体积膨胀计算量也上去了。伪近邻法FNN比较直观对每个 m 维空间里的点找它的最近邻判断这个邻域关系在把维数升到 m1 时是否仍然成立如果大量邻居不再是邻居说明 m 还不够继续加维。FNN 的问题在于需要人为设定距离阈值不同数据的阈值手感差异很大容易踩坑。Cao 方法是为了绕开阈值而设计的改进版。它构造两个指标 E1(m) 和 E2(m)E1(m) E(m1) / E(m)其中 E(m) 是 m 维空间中所有点到最近邻的距离在扩展一维后的相对变化均值。当嵌入维达到系统真实维数后E1(m) 不再明显变化就说明 m 选够了。E2(m) 是另一个指标用来区分确定性混沌和纯随机序列确定性系统的 E2 会在某些 m 处明显偏离 1而随机序列的 E2 基本保持在 1 附近。2.4 源码项目里常见的主流程预处理 → 定 τ → 定 m → 重构 → 绘图看 straightxx8 这类源码项目时不要被文件数量吓到核心调用链基本就是一条直线。首先是数据预处理去趋势、去均值、归一化这一步决定后续互信息和距离计算是否稳定。然后调互信息函数算 τ调 Cao 方法函数算 m再把原始序列按 [x(i), x(iτ), ..., x(i(m-1)τ)] 重组最后画三维轨迹或者导出点云做特征提取。源码好不好用的判断标准也就在这些环节里预处理是否考虑非平稳、参数估计是否允许人工覆盖、Cao 曲线和互信息曲线是否可视化输出。如果源码只是把函数写出来但没给你画曲线的入口大概率你要自己补一段绘图代码否则参数选择就成了玄学。3. 把 PSR 三维重构跑起来的 Python 代码从互信息定 τ 到 Cao 法定 m 再到 3D 轨迹3.1 预处理去趋势、去均值、归一化避免“趋势带飞吸引子”相空间重构对非平稳成分极度敏感。传感器信号里最常见的干扰是线性漂移和直流偏置如果不去掉重构后的轨迹会整体飘移甚至呈现向外螺旋的假象。第一个函数先做预处理import numpy as np from scipy import signal def preprocess_series(raw): # 去掉线性趋势加速度、压力、脉象等信号常带低频漂移 x signal.detrend(raw, typelinear) # 去直流重构前序列均值必须归零否则轨迹偏离原点 x x - np.mean(x) # 归一化统一量纲避免距离计算偏向幅值大的坐标轴 std np.std(x) if std 1e-12: raise ValueError(信号方差接近零无法进行相空间重构) x x / std return x逻辑说明signal.detrend默认用最小二乘拟合一条直线并扣除适合传感器标定漂移如果数据有弯曲的基线漂移要改用高通滤波或多项式去趋势。去均值是为了让重构轨迹的中心回到原点否则三维图里所有点会被整体抬升影响后续距离计算。归一化把信号变成零均值单位方差这是互信息和欧氏距离计算的统一基准。参数说明typelinear是最常用选项如果信号本身是周期平稳的也可以不归一化但距离阈值类算法伪近邻、关联维必须归一化。经验之谈做 Lorenz 或 Rossler 仿真验证时归一化可做可不做实际用振动台、心电、脑电数据时这步不做必翻车。3.2 互信息函数实现与 τ 自动选取代码预处理完成后先做延迟时间估计。下面的函数用二维直方图估计互信息再自动找第一个局部极小值def mutual_info(x, tau, bins16): # 取原序列和延迟序列长度对齐 n len(x) - tau a x[:n] b x[tau:] # 联合直方图分箱数 bins 决定概率估计的分辨率 p_ab, _, _ np.histogram2d(a, b, binsbins) p_ab p_ab / p_ab.sum() p_a p_ab.sum(axis1, keepdimsTrue) p_b p_ab.sum(axis0, keepdimsTrue) pe p_a * p_b # 只统计概率非零的格子负数被 log(0) 拦住 mask (p_ab 0) (pe 0) return float(np.sum(p_ab[mask] * np.log(p_ab[mask] / pe[mask]))) def pick_tau(x, max_tau80, bins16): taus np.arange(1, max_tau 1) mis [mutual_info(x, int(t), bins) for t in taus] # 查找第一个局部极小值比前一个小且不大于后一个 for i in range(1, len(mis) - 1): if mis[i] mis[i - 1] and mis[i] mis[i 1]: return int(taus[i]), taus, np.array(mis) # 没有明显局部极小值时退回到全局最小点 return int(taus[np.argmin(mis)]), taus, np.array(mis)逻辑说明互信息把序列值映射到二维直方图如果 τ 取 0a 和 b 完全相同互信息最大随着 τ 增加互信息下降第一个极小值说明这里的信息重复最小。用“第一个”局部极小值而不是全局最小值是因为全局最小往往出现在 τ 很大、信号已经失去相关性的地方没有物理意义。如果曲线只有一个宽谷程序会退化到全局最小值但这种情况通常暗示数据本身近似随机重构意义不大。参数说明bins建议取 16 到 64。数据长度在几千点时用 16几万点以上可以到 32 或 64分箱太小会抹平非线性结构分箱太大会让大量格子概率为 0互信息估计方差变大。max_tau我习惯设为信号主周期的三分之一到二分之一如果不知道主周期直接取 50 到 100 也够用但注意太大会让互信息计算量明显上升。3.3 Cao 方法实现与 m 判定代码τ 定完后做嵌入维估计。Cao 方法具体实现如下核心是用 KDTree 加速最近邻查找from scipy.spatial import cKDTree def cao_dimension(x, tau, max_dim12): N len(x) E np.zeros(max_dim 2) # E[d] 对应 d 维空间 Es np.zeros(max_dim 2) # 辅助指标 for d in range(1, max_dim 2): # 留出下一维的分量有效点数要减掉 d*tau n N - d * tau Y np.empty((n, d)) for col in range(d): Y[:, col] x[col * tau: col * tau n] # cKDTree 找最近邻k2 时第一个是自身第二个是最近邻 tree cKDTree(Y) dist, idx tree.query(Y, k2) neighbor idx[:, 1] d_dist dist[:, 1] # 扩展一维后的末端分量差 tail_i x[d * tau: d * tau n] tail_j x[neighbor d * tau] with np.errstate(divideignore, invalidignore): a np.abs(tail_i - tail_j) / d_dist # 最近邻距离为 0 时比值无意义按 1 处理避免均值偏爆 a[~np.isfinite(a)] 1.0 E[d] np.mean(a) Es[d] np.mean(np.abs(tail_i - tail_j)) E1 E[2:] / E[1:-1] # E1[d-1] 对应嵌入维 d E2 Es[2:] / Es[1:-1] # E2 用于区分混沌与随机 dims np.arange(1, max_dim 1) return dims, E1, E2 def pick_m(dims, E1, tol0.1): # 启发式取第一个 E1 增长趋于平缓的维度实际建议结合画图 for i in range(1, len(E1)): if abs(E1[i] - E1[i - 1]) tol * max(0.1, abs(E1[i - 1])): return int(dims[i]) return int(dims[np.argmin(np.abs(E1 - 1))])逻辑说明Cao 方法的核心思想是看“m 维最近邻关系”在升到 m1 维后是否还成立。d_dist是当前 d 维空间里每个点到最近邻的距离tail_i和tail_j是这两个点扩展一维后的新增分量两者比值的均值就是 E(d)。E1 是前后两代 E 的比值如果系统是确定性的E1 会随着 m 增加趋于 1 或一个稳定平台此时的最小 m 就是嵌入维。E2 则是用末端分量的绝对差检验随机性——如果 E2 始终接近 1数据更可能是纯随机。参数说明max_dim设 8 到 12 就够因为工程里绝大多数吸引子的嵌入维不超过 6设太大只会让 KDTree 构建时间和内存上扬。tol0.1是启发式阈值采样点少、噪声大时要放宽到 0.2并且一定结合 E1 曲线人工确认不能无脑自动选。这段代码的时间复杂度主要被 KDTree 吃掉n 在几万点时非常快几十万点建议先隔点抽样再算。3.4 三维相空间轨迹绘制与保存点云参数拿到后重构三维轨迹就简单了。下面这个函数既返回点云用于后续计算也直接画图def reconstruct_3d(x, tau, m3): # m 可以大于 3这里保留全部维度调用方自行选择画哪些列 n len(x) - (m - 1) * tau X np.empty((n, m)) for i in range(m): X[:, i] x[i * tau: i * tau n] return X def plot_psr(X, tau, m3, stride1, save_pathNone): fig plt.figure(figsize(8, 6)) ax fig.add_subplot(111, projection3d) # stride 抽稀绘制大数据量时防止线形糊成一团 ax.plot(X[::stride, 0], X[::stride, 1], X[::stride, 2], lw0.4, alpha0.8) ax.set_xlabel(fx(t)) ax.set_ylabel(fx(t{tau})) ax.set_zlabel(fx(t{2*tau})) if save_path: plt.savefig(save_path, dpi150) plt.show()逻辑说明reconstruct_3d返回的是 n×m 的矩阵每一行是一个 m 维状态点。plot_psr默认取前三个坐标画三维线如果 m 大于 3这只是低维投影要明确这一点。stride参数很实用数据量超过 10 万点时画所有线会让渲染卡死隔 2 到 5 个点抽稀就能兼顾流畅和全局形状。参数说明stride1适合几万点以内的数据点多了改 2 到 5有严重噪声时把lw调小、alpha调低能看到轨迹内部密度分布。save_path建议每次跑都存图方便后面换参数对比时回看这是省时间的小习惯。3.5 参数怎么配采样率、窗口长度、τ/m 上限的设置建议实际调参的经验值我放到一个表里新手可以直接按这个范围起步参数建议范围判断依据备注数据长度 N5000 点起步够画完整轨迹并让统计指标平稳少于 2000 点时 Cao 曲线会剧烈抖动延迟时间 τ1 ~ 80 点互信息第一极小值采样率越高可选 τ 范围越大嵌入维 m2 ~ 12E1 曲线平台起点三维重构只是 m3 的特例互信息分箱 bins16 ~ 64数据量越大可以取越大小数据用 16大数据用 32 或 64Cao 最大维 max_dim8 ~ 12预期吸引子维数 2 即可太大让 KDTree 构建变慢绘图抽稀 stride1 ~ 10画面不糊、帧率可接受画线用 stride散点可不抽经验上采样率高的机械振动信号N 在 1 万点以上、τ 在 10 到 30 之间非常常见生物学信号比如心电、脑电τ 往往更小落在 3 到 15。如果互信息曲线没有明显极小值先回头检查预处理和原始信号质量而不是硬调max_tau。4. PSR 三维重构翻车现场5 个典型坑与排查方法4.1 轨迹贴着主对角线拉成一条细线——τ 偏小现象三维图里所有点都落在一个狭长的柱体或平面附近曲线像一条拉直的弹簧几乎看不出分层和折叠。原因延迟时间 τ 选得太小x(t)、x(tτ)、x(t2τ) 之间高度相关三个坐标几乎线性相关重构后被压回一维流形。很多项目直接用自相关法取第一次过零点对非线性系统来说这个 τ 经常偏小。解决改用互信息第一极小值定 τ并且把互信息曲线画出来对照。如果数据里有强周期成分自相关法会把周期时间错当延迟时间出现“轨迹绕圈但不是吸引子”的假象。用pick_tau得到 τ 后加大到原来的 1.5 倍和 2 倍各画一张图三张对比细线结构消失、出现横向展开的那一版才是合理结果。4.2 轨迹散成一团均匀点云看不出结构——τ 偏大或噪声过大现象重构后没有连续轨迹感所有点均匀分布在三维空间像一团雾。原因分两种一是 τ 取得太大坐标之间信息完全独立动力学关联被打散二是原始数据信噪比太低随机噪声把吸引子结构整个淹没。这两种情况的修复方向完全不同所以必须明确区分。解决先用plot_psr在 τ 基础上降低一档看结构是否恢复。如果降一半 τ 就有明显结构说明是延迟时间过大如果怎么调都还是一团雾则回原始信号做带通滤波滤掉工频和随机毛刺再重新跑预处理。机械振动数据尤其要小心轴承磨损冲击噪声它们会让互信息曲线第一个极小值变得很浅看起来像随机序列。4.3 嵌入维 m 定错导致的折叠与“穿帮”怎么用 E1/E2 曲线确认现象m 取太小轨迹中出现大量虚假交叉三维图里线条自己和自己纠缠不清m 取太大曲线变得“肿胀”体积占比很大且拖尾明显。原因很直接m 小于系统真实嵌入维低维投影丢失了区分重叠状态所需的信息m 大于真实需要时把噪声维度也当成了有效状态维度。解决不要只看自动选择的 m要把 Cao 方法输出的 E1、E2 曲线画出来。E1 从某个 m 开始基本不再增长这个 m 才是合理嵌入维E2 如果显著偏离 1提醒你系统有确定性结构如果 E2 全程贴着 1那数据不具备明显的混沌特征三维重构再好看也只是在投影噪声。另外记住m 算出来是 5 以上时三维图只是前三个坐标的投影特征提取不能只依赖三维几何。4.4 轨迹围着原点转几圈后突然跳走——去趋势与去漂移没做干净现象轨迹前半段是收敛的吸引子画着画着突然飞向远方或整体呈现出向外螺旋的趋势。原因信号里残留低频漂移或非平稳均值偏移重构时漂移的方向分量被当成了真实动力学状态。这个问题在长时间采集的传感器数据里极其常见上午标定的基线到下午就变了。解决在互信息和 Cao 计算之前必须先做preprocess_series里的去趋势如果线性去趋势后仍有弯曲漂移改用scipy.signal.detrend配合低阶多项式或者直接做高通滤波截止频率取信号主频的 1/10。处理完后重新画互信息曲线通常漂移会明显改善。我最早做压电加速度数据时就栽在这里当时以为是系统有混沌特征排查了三天最后发现只是温度漂移。4.5 同一份数据换个视角就完全不像——三维投影视角检查技巧现象三维图在某个视角看到明显折叠旋转 90 度后结构却平淡无奇导致两个人对同一份数据得出相反结论。原因三维投影的观察角度会直接影响视觉判断特别是吸引子本身是平坦结构时轴测视角会把层叠结构隐藏起来。解决固定三到四个标准视角比如默认视角、沿 X 轴看、沿 Y 轴看、沿 Z 轴看每个视角都存一张图也可以用 PCA 对三维轨迹做主成分分析用前三个主成分重新投影这样能把最大方差方向拉直。最稳妥的验证方式是把已知 Lorenz 系统生成的数据用完全相同的参数跑一遍如果 Lorenz 的三维图结构正常、你的数据三维图观感却很怪说明问题在数据本身而不是代码或视角。5. 重构之后怎么验证有效用 Lorenz 系统做基准测试再往前走一步5.1 用 Lorenz 序列验证 τ 与 m 的选择到底准不准相空间重构代码写完后第一件事不是拿真实数据直接跑而是拿理论已知的系统试运行。Lorenz 方程是最方便的金标准用 RK4 生成一段 z 分量序列然后走完整条重构链路def gen_lorenz(sigma10.0, beta8/3, rho28.0, dt0.02, n6000): state np.array([1.0, 1.0, 1.0]) out np.zeros(n) for i in range(n): x, y, z state k1 np.array([sigma * (y - x), x * (rho - z) - y, x * y - beta * z]) k2 np.array([sigma * (y k1[1]*dt/2 - x - k1[0]*dt/2), (x k1[0]*dt/2) * (rho - z - k1[2]*dt/2) - (y k1[1]*dt/2), (x k1[0]*dt/2) * (y k1[1]*dt/2) - beta * (z k1[2]*dt/2)]) # 为节省篇幅这里只保留一阶近似实际用完整 RK4 state state dt * k1 out[i] z return out把gen_lorenz的结果接上preprocess_series、pick_tau、cao_dimension你会看到 Lorenz 系统互信息第一个极小值稳定在某个 τE1 曲线在 m3 或 4 后开始平缓三维轨迹呈现标准双蝶形。如果这套流程在自己写的代码里跑出来的不是这个形状优先怀疑代码里的索引边界和最近邻选择逻辑而不是怀疑参数。5.2 从重构轨迹抽取关联维数与最大 Lyapunov 指数的流程三维图只是定性的想要把重构结果量化两个最常见指标是关联维数和最大 Lyapunov 指数。关联维用 Grassberger-Procaccia 算法对重构轨迹两两计算距离统计小于半径 r 的距离对占比 C(r)然后对 log C(r) 与 log r 做线性回归斜率就是关联维。注意大数据量下两两距离矩阵会非常大常见做法是随机抽 2000 到 5000 个点计算。最大 Lyapunov 指数可以走 Wolf 算法的简化版先在重构轨迹里找每个点的最近邻跟踪一段时间内两个轨迹的分离速率最后对 log 距离差做线性拟合斜率就是最大 Lyapunov 指数。这一步需要调节跟踪窗长一般取重构后轨迹总长度的 1/10太大会跨过折叠区太小拟合噪声大。5.3 重构轨迹存成三维点云之后还能做什么重构后的相空间点云本身就是最直接的特征它不需要标签可以直接拿去做聚类、异常检测和分类。比如旋转机械的故障诊断把振动信号重构到三维相空间再计算点云的体积、密度分布、关联维数故障前后的差异往往比时域特征明显得多。心电和脑电分析里也可以对重构点云做可视化筛查人工确认异常节律对应的轨迹形态。最后一个习惯是我自己的血泪经验每次调参必须把互信息曲线、Cao 曲线、重构三维图三个一起存下来文件名里带上数据编号、τ 和 m这样回头排查才不用重新跑一遍。相空间重构的每一步都有取舍参数选择不是一次到位而是用已知系统验证、用可视化确认、再用特征指标定量的三级流程。希望这些经验能帮你在做 PSR 三维重构时少走几趟弯路。本文还有配套的精品资源点击获取