★ FEATURED ARTICLE 简介这份MATLAB源码包面向信号处理、非线性动力学与时间序列分析研究者围绕三维相空间重构PSR提供一套完整可运行的算法实现适用于科研复现、教学演示与工程项目开发。压缩包共11个文件以4个M脚本为核心涵盖互信息法延迟估计、FNN嵌入维计算、Lorenz混沌系统时间序列生成与三维可视化另含3张结果图片、1份txt数据、1个Markdown说明文档与License授权文件包体仅205KB。目前已有183人学习下载。通过这套代码可直观理解Takens延时嵌入定理的实操流程从原始一维时间序列到重构相空间再到吸引子形态绘制完整呈现相空间重构的关键步骤代码结构清晰参数选择部分便于替换到真实信号中进一步探索是学习混沌分析的良好参考实现。1. 一份 PSR 三维重构源码能帮你看到什么“所有代码_psr_三维重构_相空间_相空间重构_straightxx8_源码”是典型的资源站下载包一堆关键词拼起来的压缩包没有说明书也没有版本号。它的技术主线很聚焦——用相空间重构Phase Space ReconstructionPSR把一维时间序列映射到三维空间里把混沌信号藏在时间轴里的吸引子结构“画”出来。这类源码在振动故障诊断、生理信号分析、非线性时间序列预测里出现频率很高适合手里有一段实测数据、想判断它到底是随机噪声还是确定性结构或者想给分类模型造一个更好分特征的人。它解决的核心问题可以用一句话概括同一段波形时域里看不出规律升到三维空间后规则结构立刻现形。下面我按复现这类源码包的顺序把原理、算法、实现、踩坑和定量分析一次讲清。2. 相空间重构原理与参数选型τ 和 m 为什么决定三维图长什么样相空间重构在大部分人听来像玄学核心其实是一句话一段标量时间序列里藏着系统全部状态的演化轨迹。决定三维重构图好不好看的只有两个参数——延迟时间 τ 和嵌入维数 m。源码包里几乎所有子程序都在围着这两个参数转读懂了它俩任何 PSR 源码都不会再看晕。2.1 Takens 嵌入定理从一维序列恢复吸引子拓扑Takens 在 1981 年证明的嵌入定理是这套方法的根基。假设原始动力系统是 d 维的我们能观测到的只是其中一个坐标的采样序列 x(t)。构造延迟向量X(t) [x(t), x(tτ), ..., x(t(m-1)τ)]当嵌入维数 m ≥ 2d1 时重构后的轨迹与原始吸引子是拓扑等价的。翻译成人话虽然每个时刻只观测到一个数值但把“现在”和“未来几个时刻”拼成一个向量足够还原系统内部状态的演化关系。所谓三维重构就是取 m3 的特例三个坐标轴分别是 x(t)、x(tτ)、x(t2τ)。这里有个特别容易被误解的点重构坐标没有物理单位也不代表原系统里的某个物理量它只是延迟副本构成的抽象空间。所以别给坐标轴硬标“电压”“位移”之类的量。拓扑等价的意义在于几何不变量可以保留——关联维数、Lyapunov 指数这些反映系统本质的量在重构空间里算和在原系统里算结果一致这是后面做定量分析的前提。工程上d 一般未知所以 m 通常从 2 试到 10 左右看结构和指标是否稳定。如果只是想“画个三维吸引子看看”m 固定为 3 就够了。源码包里大量出现 m3 不是偷懒是可视化场景下的合理选择。2.2 延迟时间 τ 的两种算法自相关法与互信息法τ 选小三个坐标高度相关轨迹挤成一条线τ 选大三个坐标近似独立轨迹变成随机点云。自相关法和互信息法是最常用到的两种选法。自相关法算的是 x(t) 和 x(tτ) 之间的线性相关系数随 τ 的衰减常见准则取第一次降到 1/e 的位置作为 τ。优点是快缺点是只捕捉线性依赖对非线性结构不敏感。import numpy as np def autocorr_tau(signal, stop1.0 / np.e): x signal - signal.mean() n len(x) # 补零到 2n用 FFT 算线性自相关避免逐点循环 fft_x np.fft.fft(x, n2 * n) acov np.fft.ifft(fft_x * np.conj(fft_x)).real[:n] / n acov acov / acov[0] # 归一化到 lag0 时相关系数为 1 for tau in range(1, n): if acov[tau] stop: return tau return n - 1逻辑说明先减均值消除直流分量再补零做 FFT 计算自相关比逐点双重循环快几个数量级。除以 acov[0] 完成归一化阈值就直接用 1/e。对 Lorenz 这类信号在 dt0.02 时算出的 τ 通常在 5 到 15 之间和文献里“延迟时间取自相关第一次过零点附近偏小一点”的经验吻合。保守的写法是取第一次过零点但那给出来的 τ 往往偏大轨迹会明显变稀疏。自相关法的局限在于它只度量线性相关性。互信息法则能捕捉非线性依赖它把信号值域分成若干格子统计滞后 τ 的两个变量共享多少信息量取第一极小点作为 τ。def mutual_information(signal, tau, bins32): x signal[:-tau] y signal[tau:] lo, hi np.min(signal), np.max(signal) # 联合直方图固定使用全序列值域保证不同 tau 之间可比 cxy, _, _ np.histogram2d(x, y, binsbins, range[[lo, hi], [lo, hi]]) n cxy.sum() pxy cxy / n px pxy.sum(axis1) py pxy.sum(axis0) mi_val 0.0 for i in range(bins): for j in range(bins): if pxy[i, j] 0: mi_val pxy[i, j] * np.log(pxy[i, j] / (px[i] * py[j])) return mi_val def mi_first_min(signal, tau_max80, bins32): vals [mutual_information(signal, t, binsbins) for t in range(1, tau_max 1)] for i in range(1, len(vals) - 1): if vals[i] vals[i - 1] and vals[i] vals[i 1]: return i 1 # 索引 i 对应 tau i1 return tau_max参数说明bins 取 32 是常见折中数据总量少于几千点时降到 16否则联合直方图大量格子为零互信息抖动很厉害。tau_max 要覆盖信号的一个主周期dt0.02 的 Lorenz 轨道时间常数在 1 秒量级tau_max 取 80 足够。代码里返回的是第一个局部极小点不是全局最小点这是 Fraser-Swinney 方法的经典约定。两套算法结果不一致时怎么办比如自相关给 8、互信息给 15先画互信息曲线看第一极小是否明显再在两值之间取偏大者做可视化。稍大的 τ 能把轨迹拉开、看到更多折叠结构如果差异超过 3 倍多半是信号有趋势或周期性太强先去趋势再说。2.3 嵌入维数 m 的确定从伪近邻到“够用就好”如果只是三维可视化这部分可以跳过。但源码包通常还带 G-P 算法或伪近邻法说明作者意图不止画图。伪近邻的思路在 m 维空间里一个点的大部分近邻应该是“真邻居”如果升到 m1 维后原本的近邻跑远了说明那些是低维投影造成的假邻居。m 从 1 递增伪近邻比例降到接近 0 时的 m 就是合适嵌入维。G-P 算法从另一个方向逼近在重构空间里统计距离小于 r 的点对比例得到关联积分 C(r)log-log 坐标下无标度区的斜率就是关联维数 D2。随着 m 增大确定性混沌系统的 D2 会饱和在某个值附近如果 D2 一直涨信号大概率是随机噪声。仅这一条就常被用来区分“混沌”和“纯随机”。工程选型建议画三维图用 m3估算关联维数或 Lyapunov 指数用 m5 到 7 起步高维系统通常要 m≥8。m 不是越大越好——样本量固定时空间维数越高数据越稀薄距离估计全部失真。经验上要求重构后的点数 n-(m-1)τ 至少是 m 的 10 倍否则后面算关联维数基本是噪声。3. 从源码包到最小复现Lorenz 信号的三维相空间重构拿到这类源码包最常见的做法是先把环境配干净用一段已知答案的混沌信号把流程跑通再换自己的数据。不要一上来就上真实信号因为真实信号里的噪声和趋势会让“图不对”时无法判断是自己错了还是数据本身有问题。3.1 造一段已知答案的测试信号Lorenz 系统from scipy.integrate import solve_ivp def lorenz(t, state, sigma10.0, rho28.0, beta8.0 / 3.0): x, y, z state return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z] dt 0.02 t_end 120 t_eval np.arange(0, t_end, dt) sol solve_ivp(lorenz, [0, t_end], [1.0, 1.0, 1.0], t_evalt_eval, methodRK45, rtol1e-8) x sol.y[0] x x[2000:] # 丢掉前 40 秒瞬态 print(f剩余点数: {len(x)})逻辑说明Lorenz 方程在 sigma10、rho28、beta8/3 的经典参数下处于蝴蝶混沌区初值随便给只要不落在平衡点附近就行。积分完成后把前 2000 点丢弃因为从初值飞到吸引子上的过渡段会在重构图里多出一条“飞线”。两个参数要记住dt 是采样间隔直接决定 τ 的物理含义rtol1e-8 防止数值误差让轨迹跳到另一个分支。真实数据没有积分这一步但一定有采样率建议一开始就把 τ 的离散值换算成物理时间。提示真实信号做相空间重构前先确认采样率和主频带。采样率过高时先降采样否则重构点数暴涨图也卡τ 的物理意义也容易算错。3.2 相空间重构核心实现与三维可视化def psr_reconstruct(signal, tau, m3): n len(signal) rows n - (m - 1) * tau if rows 0: raise ValueError(n-(m-1)*tau 为负数据太短或参数太大) mat np.empty((rows, m)) for i in range(m): mat[:, i] signal[i * tau : i * tau rows] return mat tau 12 mat psr_reconstruct(x, tau, m3) print(mat.shape) # (rows, 3) import matplotlib.pyplot as plt fig plt.figure(figsize(8, 6)) ax fig.add_subplot(111, projection3d) ax.plot(mat[:, 0], mat[:, 1], mat[:, 2], lw0.5, colorsteelblue) # 三个轴按实际数据范围等比防止图形被压扁 ax.set_box_aspect((np.ptp(mat[:, 0]), np.ptp(mat[:, 1]), np.ptp(mat[:, 2]))) ax.set_xlabel(x(t)) ax.set_ylabel(x(tτ)) ax.set_zlabel(x(t2τ)) ax.view_init(elev20, azim45) plt.show()逻辑说明psr_reconstruct 返回 rows×3 矩阵第 0 列是原序列第 1 列滞后 12 个采样点第 2 列滞后 24 个。等价于从第 0 个原始点开始以 τ 为步长取三个元素构成第一个三维向量然后逐点滑动。画图用 plot 而不是 scatter几千个点只有在连线模式下才能看到连续的折叠结构线宽 0.5 避免蝶翼两侧互相糊成一片。set_box_aspect 是三维图不被压扁的关键很多流传的源码包里没有这一句蝴蝶会被硬拉成飞饼。view_init 固定视角后面做参数对比时才不会换一个角度就换了一张图。mat 行数超过两万时先隔点抽样再画mat[::2] 丢一半点速度翻倍且视觉几乎不变。这也是源码包里经常出现的处理不是偷数据是控制渲染量。3.3 把 τ 的自动估计接进主流程tau_corr autocorr_tau(x, stop1.0 / np.e) tau_mi mi_first_min(x, tau_max80, bins32) print(f自相关法 tau{tau_corr}, 互信息法 tau{tau_mi}) tau tau_mi if tau_mi is not None else tau_corr mat psr_reconstruct(x, tautau, m3) fig.suptitle(fLorenz, tau{tau}, m3, dt0.02)参数说明自相关和互信息结果不一致时我一般先看一眼互信息曲线确认第一极小点旁边没有毛刺再决定是否改用 tau_corr。自动估计的 τ 只配当起点不配当标准答案——用下一章的参数扫描验证过才算数。4. 相空间重构常见问题排查五个翻车现场的现象、原因与对策相空间重构的坑都很隐蔽因为程序不会报错“τ 选错了”。下面五条按出现频率排序每一条都值得在自己数据上对照一遍。4.1 现象重构轨迹全部挤在空间对角线附近三维图是一条细长的对角线或者紧紧贴在一个平面上看不到蝴蝶的折叠。这是最典型的翻车现场。原因有二τ 太小三个坐标分量数值几乎相等或信号未去均值、带趋势趋势项把轨迹拉成一条斜线。经验法则凡是吸引子看起来像个棒子先怀疑 τ再怀疑预处理。解决先做预处理再去调 τ。from scipy.signal import detrend x_clean detrend(x - x.mean()) tau_new mi_first_min(x_clean, tau_max80, bins32) mat psr_reconstruct(x_clean, tau_new, m3)逻辑说明detrend 默认去掉线性趋势去均值消掉直流分量。对缓慢漂移的实测信号这两步有时比调 τ 更关键。处理完再跑互信息法τ 往往会变大一点轨迹也会从对角线上“松开”。4.2 现象改变视角后吸引子结构完全变样同一份数据elev20 时看是蝴蝶elev70 时看成一团乱线两个人截图对比得出的结论完全相反。原因三维图本质是二维投影视角和坐标缩放都会扭曲视觉结构。尤其缺了 set_box_aspect 时三个轴按各自范围独立拉伸真实几何比例被破坏。解决固定视角加等比盒子。检查绘图代码里有没有 set_box_aspect 和 view_init 两行没有就补上。所有参数对比统一用同一视角保存图片时把视角参数写进文件名否则截图无法追溯。这是血泪经验看吸引子形状必须先固定视角否则等于看图猜谜。4.3 现象数据截断后吸引子结构剧变用前一半数据画图是一个环用后一半画是另一个环掐头去尾再看形状大变。原因数据里混入了瞬态段或者系统状态本身发生了迁移。Lorenz 测试信号里常见的是初值飞线真实传感器数据里常见的是缓慢漂移造成的状态切换。解决先定位瞬态段丢掉再用滑动窗口截取稳态段。粗略判断稳态的办法是计算每 200 点窗口的质心质心在三维空间里的偏移超过坐标范围的 10%就得重新选段。def check_stationary(mat, win200, ratio0.1): center mat.mean(axis0) spans np.ptp(mat, axis0) for start in range(0, len(mat) - win, win): seg_center mat[start:startwin].mean(axis0) if np.any(np.abs(seg_center - center) / spans ratio): return False return True逻辑说明质心漂移是吸引子结构不稳的直接信号。返回 False 时别急着调 τ先换数据段。这个函数对真实信号尤其有用它能直接指出哪一段不属于同一个动力学状态。4.4 现象τ 选太大轨迹变成稀疏点云三维图是一堆悬浮的散点看不出连续轨道像噪声而非吸引子。原因互信息法自动选 τ 时取错了极小点常见的是第一极小不明显、代码误取第二极小或者信号周期性太强自相关法的 1/e 准则直接失效。解决把互信息曲线画出来人工确认第一个极小点。import matplotlib.pyplot as plt taus np.arange(1, 80) mis [mutual_information(x, t, bins32) for t in taus] plt.plot(taus, mis) for i in range(1, len(mis) - 1): if mis[i] mis[i - 1] and mis[i] mis[i 1]: print(局部极小 tau , i 1) plt.show()逻辑说明互信息函数单个 τ 的复杂度是 O(bins²)80 个 τ 跑下来也就几十毫秒放心循环。看到曲线上低于均值的第一处凹陷那个位置才是合理 τ不是整条曲线的最低点。如果曲线第一个极小出现在 tau1说明数据可能本身采样过密或周期性过强先降采样再重构。4.5 现象两次运行结果的坐标范围不一致无法对比昨天画的吸引子范围是 [-20, 20]今天变成 [-15, 15]形状看着也不一样但代码一行没改。原因数据段起点变了、去趋势的位置变了、τ 变了图上却看不出参数差异。这不是算法错误是复现管理问题。解决每次重构输出时记录数据段起止索引、τ、m、坐标范围。常见做法是存一个 JSON或者直接编进文件名。具体模板放在最后一章这里先记住结论没有参数快照的重构结果等于没有刻度尺的图纸。5. 三维相空间重构的下游定量分析从看图到算数三维图只能让你“看着像”要说服别人、要落到项目里得把“像蝴蝶”变成“D2≈2.05”这种可复现的数值。这章讲最常用的两步。5.1 关联维数G-P 算法把吸引子形状变成一条饱和曲线from scipy.spatial.distance import pdist def correlation_integral(mat, r): n mat.shape[0] if n 8000: idx np.random.choice(n, 8000, replaceFalse) mat mat[idx] n 8000 dists pdist(mat, metriceuclidean) pairs np.sum(dists r) return 2.0 * pairs / (n * (n - 1))逻辑说明pdist 的复杂度是 O(n²)几万点会直接吃爆内存所以超过 8000 行先随机抽样。这里抽的是重构轨迹的行也就是相空间里的点不影响几何结构只降低精度。r 的扫描用对数等分rs np.geomspace(0.01, 50, 40) mat3 psr_reconstruct(x, tau, m3) cs np.array([correlation_integral(mat3, r) for r in rs]) # 无标度区经验范围C(r) 在 0.01 到 0.5 之间 mask (cs 0.01) (cs 0.5) d2 np.polyfit(np.log(rs[mask]), np.log(cs[mask]), 1)[0] print(fD2 ≈ {d2:.3f})参数说明mask 选的是 C(r) 在 0.01 到 0.5 之间的点太小的 r 区域是离散点噪声太大则进入饱和段。Lorenz 的 D2 文献值约 2.05算出来在 1.9 到 2.2 之间都算正常。如果差得远不要怀疑算法回去查 τ 和数据长度——这是祖传的调参顺序。5.2 用重构轨迹做状态识别的两个特征工程落地时很多人不关心 D2只想要一个能区分“正常”和“异常”的特征。三维重构轨迹可以抽出几个比时域统计量更敏感的特征。def psr_features(mat): cov np.cov(mat.T) eig np.linalg.eigvalsh(cov) var_ratio np.max(eig) / np.sum(eig) # 主方向方差占比 seg np.diff(mat, axis0) arc_len np.sum(np.linalg.norm(seg, axis1)) # 轨迹总弧长 volume np.prod(np.ptp(mat, axis0)) # 轨迹占据的空间体积 return var_ratio, arc_len, volume逻辑说明var_ratio 反映轨迹在三维空间里铺得广不广结构越扁此值越高arc_len 是轨道在吸引子上绕的总长度数据段相同长度时反映绕圈密度volume 是三个轴范围的乘积粗估吸引子占据空间大小。这三个量对状态切换比均值方差敏感得多。常见做法正常工况取一段数据算一组特征异常工况取另一段算一组喂给阈值判断或 SVM。但要注意边界特征对数据长度和预处理极其敏感对比时必须用相同的数据段长度和相同的 τ。比如旋转机械的振动信号转速一变特征整体漂移得先按转速分段再对每段单独重构。5.3 参数扫描τ 从 1 到 30m 从 3 到 6哪个组合最稳看单张三维图选 τ 还是容易犯主观更可靠的办法是跑参数扫描看 D2 对参数的稳定性。results [] for m in [3, 4, 5, 6]: for tau in range(1, 31): mat_t psr_reconstruct(x, tau, mm) rs_t np.geomspace(0.01, 50, 30) cs_t np.array([correlation_integral(mat_t, r) for r in rs_t]) mask_t (cs_t 0.01) (cs_t 0.5) if mask_t.sum() 3: continue d2_t np.polyfit(np.log(rs_t[mask_t]), np.log(cs_t[mask_t]), 1)[0] results.append((m, tau, d2_t))参数说明这组循环是 4×30120 次 G-P 计算每次抽样 8000 点普通笔记本几分钟内能跑完。选出 D2 随 m 饱和、且对 τ 变化不敏感的区域那个 τ 就是稳定工作点。“对 τ 不敏感”本身就是重要信号——如果 D2 随 τ 剧烈抖动说明数据长度不足或系统根本不是单个吸引子继续调参数没有意义。注意无标度区的 mask 范围0.01~0.5只在数据量足够时有效。数据少于 1000 点时不要强行算 D2结果没有统计意义。6. 给重构结果留个状态快照文件名就是后悔药6.1 参数快照模板与自解释命名写完图或算出 D2 后第一件事是把参数固化下来。τ12、m3 这个组合到底对应哪段数据、采样间隔多少、视角多少度没有这些三维图只是张无法复现的插图。meta { source: lorenz_x, start_idx: 2000, end_idx: 6000, dt: 0.02, tau: 12, m: 3, elev: 20, azim: 45, range: [float(mat.min()), float(mat.max())], } import json with open(recon_meta.json, w) as f: json.dump(meta, f, indent2)参数说明range 记录三个轴合并后的最小最大值再次绘图时用它统一坐标范围。文件名用“tau12_m3_i2000_6000.png”这种自解释命名比“重构结果.png”强得多。JSON 里再存一份完整参数图丢了还能重建。6.2 换数据前的内置校验我被这类问题坑过不止一次同一份振动数据上午下午各跑一遍画出的图一个宽一个扁最后发现只是一个 τ 用 8、一个用 10还没人记得谁用了哪个。从那以后所有重构实验一律带参数快照。一个实用的验证习惯把代码换到陌生数据上之前先在 Lorenz 上复现 D2≈2.05确认整个代码通道没问题再碰真实数据。真实数据算出的 D2 落在 1.1 到 2.9 之间通常说明有确定性结构接近整数或半整数更有说服力如果 D2 大于 4 或找不到无标度区先怀疑数据而不是算法。真正常规、能反复用、能对比的相空间重构流程一定长着“参数看得见、视角固定、坐标等比”的样子。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
