简介面向混沌时间序列分析研究者与MATLAB使用者提供基于Cao方法求取嵌入维数与时间延迟的完整实现聚焦Rossler系统的相空间重构与互信息计算。压缩包共9个文件以6个m脚本为核心覆盖数据生成、互信息计算、Cao算法主流程和相空间重构两套计算流程便于对照验证另附2个txt中间结果和1个fig图形文件用于核对输出整包仅8KB运行轻量、易于二次修改。已有356人学习下载。Cao方法通过互信息曲线极小值确定最佳延迟并结合嵌入维数判断准则完成状态空间重构脚本保留从Rossler时间序列生成、延迟量选取到维数验证的完整链路结果图可帮助理解参数变化对重构质量的影响。替换输入序列后可迁移至其他混沌动力学系统的延迟与维数估计适合科研入门、课程实验或复现经典混沌分析案例。1. 嵌入维数这道坎用 Cao 方法解 Rossler 系统的重构参数从一段混沌时间序列里恢复动力学信息第一步就卡在相空间重构上嵌入维数 m 和时间延迟 τ 取多少直接决定后面算 Lyapunov 指数、关联维数、做预测的结果。我用自相关法取过延迟也被 FNN 方法的阈值折腾过——阈值一改维数就变基本靠玄学。后来换成 Cao 方法配合互信息找延迟才把这组参数稳定住。这份 cao-method MATLAB 资源就是解决这个问题的它自带 Rossler 系统生成、相空间重构、Cao 维数扫描脚本输出 m1/m2 值适合做混沌动力学、非线性时间序列和信号重构方向的从业者。你不需要自己从头推导改改参数就能复现。2. 原理先立住互信息定延迟、Cao 法定维数为什么能避开阈值2.1 时间延迟不是越小越好互信息找第一个极小值相空间重构是把一维标量序列 x(t) 变成一组 m 维状态向量 [x(t), x(tτ), ..., x(t(m-1)τ)]目的是还原原系统的吸引子拓扑。这里 τ 的选择非常微妙τ 太小相邻坐标几乎一样吸引子被压扁在对角线附近重构出来的点都在一条线上τ 太大相邻坐标之间已经失去关联变成随机散射动力学信息被噪声盖住。自相关函数 R(τ) 是很多人第一个会试的工具通常取它第一次过零点的位置当延迟。问题在于自相关只刻画线性相关性而 Rossler 系统的动力学本质是非线性的线性无关不等于统计独立。我实际跑下来的体验是自相关给出的 τ 往往偏小重构出的 Rossler 吸引子明显挤压横幅比正常要扁。互信息是更适合非线性系统的工具。它的定义是I(τ) Σ p(x(t), x(tτ)) · log[p(x(t), x(tτ)) / (p(x(t)) · p(x(tτ)))]它衡量的是两个变量之间的统计依赖程度不限于线性关系。p(x(t)) 和 p(x(tτ)) 是一维概率p(x(t), x(tτ)) 是联合概率。实际计算时先把数据范围分成若干等宽箱子典型取 16 或 32统计数据点落在每个箱的频率再用这些频率近似概率。取 I(τ) 第一次出现局部极小值的 τ作为重构延迟。这个位置意味着在保留动态依赖的前提下坐标间的冗余信息最少。代码层面我一般把最大搜索范围设为总长度的 5%10%比如 20000 点就搜到 200 左右。箱子数不是越多越好数据量只有几千点时设 32 个箱子会让很多格子出现零频率估计出的互信息会抖动得很厉害。2.2 Cao 方法用 E1/E2 比值替代 FNN 阈值确定延迟 τ 之后剩下的嵌入维数 m 是第二个坎。最常用的 FNN假近邻法思路是在 m 维空间找每个点的最近邻然后看这两个点在 m1 维空间里是否还是近邻如果不再是说明原来那个邻居是因为维数不够而造成的“假邻居”。问题在于判断“还是不是近邻”需要给定一个相对阈值 R_tol 和绝对阈值 A_tol这两个阈值完全是主观的。阈值设 10 或设 15同一个序列能算出两个不同的维数。Cao 方法Cao, 1997就是为了消掉这个主观阈值。它的做法是这样对每个候选维数 d先把序列重构到 d 维和 d1 维然后在 d 维空间里找每个状态点的最近邻要排除自己接着计算这两个邻近点在 d1 维空间里的距离与它们在 d 维空间里距离的比值a(i,d) || X_i(d1) - X_j(d1) || / || X_i(d) - X_j(d) ||其中 j 是 i 的最近邻索引。对所有状态点取平均得到 E1(d)E1(d) (1/(N - dτ)) · Σ a(i,d)当 d 小于真实嵌入维数时维度每增加一维原来被压扁的近邻会被拉开E1 会明显变化当 d 已经达到真实维数之后再增加维度最近邻的相对距离不再系统性变大E1 进入平台。更常用的判据是 E2E2(d) E1(d1) / E1(d)当 E2 从明显大于 1 的状态掉到 1 附近并保持稳定就说明已经到达真实嵌入维数附近。Cao 方法不需要设阈值只看曲线趋势这就是它在工程里比 FNN 顺手的原因。需要注意Cao 方法输出的是“最大需要的维数”而不是“系统真实维数”。Rossler 系统的吸引子分形维数大概在 2.1 左右但 Takens 嵌入定理要求 m ≥ 2D1所以理论下限是 5 到 6实际重构时因为数据有限Cao 给出的 m 通常在 36 之间。你最后用的是 m 而不是 D这一点别混。2.3 为什么选 Rossler 当测试对象以及 m 与 τ 的连带关系Rossler 系统是经典三维混沌模型方程是 dx/dt -y-zdy/dt xa·ydz/dt bz·(x-c)取 a0.2、b0.2、c5.7 时处在混沌状态。它比 Lorenz 系统简单吸引子形态清楚也更容易目测验证所以最适合当 Cao 方法的测试对象。这份资源里所有脚本都围绕这条主线设计你可以先用它跑通流程、看懂输出再把自己的数据替换到 data 变量里。还有一点容易被忽略延迟 τ 和嵌入维数 m 不是两个独立变量。τ 取小的时候Cao 的 E2 曲线往往更早进入平台给人“维数很小”的错觉τ 取太大数据点之间相关性太弱E1 曲线上会出现不必要的波动。所以我的习惯是先算互信息把 τ 固定下来再跑 Cao如果算出 m 之后发现相图挤压或发散再回头微调 τ 扫一次看 m 是否变化。这不是严谨的理论结论但实际处理混沌数据时很实用。3. 拆 MATLAB 源码data_Rossler、reconstitution、cao_Single 在干什么3.1 data_Rossler.m先解决“数据从哪来”资源包里的 data_Rossler.m 是入口脚本负责生成用于验证的 Rossler 时间序列。直接用 MATLAB 的 ode45 积分即可关键是把输出处理成干净的标量时间序列。% data_Rossler.m % 生成 Rossler 系统的 x 分量时间序列 clear; clc; a 0.2; b 0.2; c 5.7; % Rossler 混沌参数 dt 0.01; % 采样步长 T 0:dt:200; % 总仿真时长 200 个时间单位 x0 [0; 0; 1]; % 初始状态 % 用 ode45 积分 Rossler 方程 [t, x] ode45((t,x) [ -x(2)-x(3); x(1)a*x(2); bx(3)*(x(1)-c) ], T, x0); % 去掉前 5000 点瞬态取 x 分量作为观测序列 data x(5000:end, 1); save(rossler_data.mat, data, dt);这段代码的逻辑不复杂先定义参数再用 ode45 积分最后保存数据。需要注意两个细节。第一c5.7 是 Rossler 系统进入混沌的经典取值不少人把 c 改成 4 发现系统进入周期窗口然后误以为 Cao 方法算错了维数。第二前 5000 点必须丢弃因为初值 [0;0;1] 到吸引子上的轨迹有一段瞬态过程保留它会让互信息和 E1 曲线都带上“污染段”。这里的 dt 不是数值积分的内部步长ode45 是自适应步长dt 只决定输出点的间隔。对 Rossler 系统dt0.01 已经足够密每个振荡周期大约能采 60100 个点。如果你把 dt 改到 0.05互信息曲线会明显变粗糙有可能错过第一个极小值点。采样总时长 200 个时间单位对应 20000 个数据点去瞬态后剩 15001 点这个量对 Cao 方法来说是合格的。3.2 reconstitution.m相空间重构核心函数reconstitution.m 是整个包的地基Cao 方法的每一步都要调用它。它的输入是原始序列、嵌入维数 m、延迟 tau输出是重构后的相空间矩阵。function X reconstitution(data, m, tau) % data: 原始标量序列列向量 % m: 嵌入维数 % tau: 延迟时间 % X: 重构后的相空间矩阵每一行是一个状态点 N length(data); num_points N - (m - 1) * tau; % 能构造出的状态点数量 X zeros(num_points, m); for i 1:num_points for k 1:m X(i, k) data(i (k - 1) * tau); end end end逻辑很简单num_points 是最多能取到的状态点数因为最后一个状态点要取到 data(N)所以索引是 i(m-1)τ ≤ N。循环里 i 表示状态点序号k 表示当前点在这一行的第几列。输出矩阵的行数是 num_points列数是 m这也是后续 knnsearch 能直接使用的格式。我见过不少人在这个函数上报错典型原因是把 num_points 和 m 写反。这里再强调一次XO reconstitution(data, m, tau) 之后size(X) 第一维是 N-(m-1)τ 而不是 m。后面做 Cao 计算时如果发现索引越界或者 mean(a) 结果异常请先回头打印一下 size(X)。3.3 cao_Single.m 与 doubleCao.m从单次扫描到全维曲线cao_Single.m 是核心计算函数它对每个候选维数 d 计算 E1 和 E2。我按资源包的名字和输出文件推断作者把“doubleCao”理解成同时计算 E1、E2 两条曲线的完整流程。这里给出一个能直接跑的兼容版本function [E1, E2] cao_Single(data, tau, maxDim) % 对每个嵌入维数 d 计算 E1(d) 和 E2(d) N length(data); E1 zeros(1, maxDim); for d 1:maxDim Xd reconstitution(data, d, tau); % d 维相空间 Xd1 reconstitution(data, d 1, tau); % d1 维相空间 % 找每个点的最近邻K2 时第 1 列是自己第 2 列才是真邻居 idx knnsearch(Xd, Xd, K, 2); sum_a 0; numP size(Xd, 1); for i 1:numP j idx(i, 2); dist_d norm(Xd(i, :) - Xd(j, :)); dist_d1 norm(Xd1(i, :) - Xd1(j, :)); sum_a sum_a dist_d1 / (dist_d eps); end E1(d) sum_a / numP; end E2(1) 1; for d 2:maxDim E2(d) E1(d) / E1(d - 1); end end逻辑说明对每个 d先在 d 维空间找最近邻 j再把这个邻居索引套到 d1 维空间里算距离比。这里的关键是 K2因为 K1 时最近邻就是自己距离比为 1E1 曲线会变成一根平线看起来像已经收敛实际是假象。注意 E2 的定义是 E1(d1)/E1(d)我在代码里用后一个 E1 除以前一个 E1 来填充 E2(d)方向是一致的。判断规则是当 E2 从明显大于 1 的位置掉下来然后在 1 附近持续波动这个转折点对应的 d 就是 Cao 方法建议的嵌入维数。doubleCao.m 是调用 cao_Single 并输出结果的上层脚本。原包里应该包含类似下面的逻辑只是文件命名和保存格式可能有差异% doubleCao.m 求 Rossler 系统嵌入维数 clear; clc; load(rossler_data.mat); tau 8; % 由互信息曲线第一个极小值得到 maxDim 15; [E1, E2] cao_Single(data, tau, maxDim); % 找 E1 进入平台的起始维数常见判据变化幅度小于 0.05 m1 find(abs(diff(E1)) 0.05, 1); % 找 E2 首次稳定在 1 附近的维数允许 0.2 的容差 m2 find(abs(E2 - 1) 0.2, 1); fprintf(m1 %d, m2 %d\n, m1, m2); % 保存结果 save(m1m2的值.txt, m1, m2, -ascii); save(m1的值.txt, m1, -ascii); % 画图并保存 figure; plot(1:maxDim, E1, o-, 1:maxDim, E2, s-); legend(E1, E2); xlabel(嵌入维数 d); ylabel(E1 / E2); savefig(m1.fig);m1 对应 E1 曲线进入平台的起点m2 对应 E2 稳定在 1 附近的起点。两者通常一致如果差超过 1优先信 m2因为 E1 平台的“判稳”比 E2 更依赖容差设置。你在资源包里看到的 m1.fig 就是这个图m1m2 的值.txt 是这两个数字。3.4 包里没有互信息函数时自己补一份文件列表里没有专门的互信息计算脚本但标题里写了“延迟 互信息”说明作者默认这一步已经单独做过了。我在这里补一个直方图分箱的简化实现方便你把流程串起来。function [tau, mi] mutual_info_delay(data, maxTau, bins) % 用等宽直方图分箱估计互信息返回曲线和第一个极小值 data data(:); N length(data); edges linspace(min(data), max(data), bins 1); dataIdx discretize(data, edges); numObs N - maxTau; mi zeros(1, maxTau); for tau 1:maxTau pxy zeros(bins, bins); for i 1:numObs pxy(dataIdx(i), dataIdx(i tau)) ... pxy(dataIdx(i), dataIdx(i tau)) 1; end pxy pxy / numObs; px sum(pxy, 2); % bins x 1 py sum(pxy, 1); % 1 x bins logTerm log(pxy ./ (px * py) eps); mi(tau) sum(pxy(:) .* logTerm(:)); end % 找第一个局部极小值数据质量差时可直接用全局最小 [~, tau] min(mi); end逻辑说明先对数据分箱再用二维频数近似联合概率 pxy一维频数近似边缘概率 px 和 py。px * py 是外积得到 bins×bins 的独立假设概率矩阵两者比值的对数就是互信息贡献。最后相加得到 I(τ)。参数上bins 通常取 16如果你的数据只有 3000 点用 8 到 10 个箱子更稳。maxTau 我一般取 N/20 左右太大会让 numObs 太小尾部统计不可靠。代码里取的是全局最小值严格做法应该找第一个局部极小点两者在 Rossler 这类清晰系统上差别不大。4. 实操全流程从生成 Rossler 序列到读出 m1/m2 值4.1 搭建目录并确认数据量拿到资源包后先把所有 .m 文件放进同一个目录MATLAB 里 cd 到该目录。然后运行 data_Rossler.m观察工作区里 data 的长度。这里有一个检查数据的习惯不要急着跑后续脚本先确认 data 是列向量、长度在 10000 以上。如果长度只有几百后面互信息和 Cao 都会给出离谱结果。cd(你的路径/cao-method); run(data_Rossler.m); whos data dt length(data) % 应该输出 15001运行后 rossler_data.mat 会出现在当前目录。这个文件就是后续所有计算的唯一数据源。如果你要处理自己的数据直接在这一步之后把 data 替换掉并把 dt 改成你自己的采样间隔。4.2 用互信息确定延迟 tau按照 3.4 的函数在命令行执行[tau, mi] mutual_info_delay(data, 200, 16); figure; plot(1:200, mi); xlabel(tau (采样点数)); ylabel(I(tau)); title(Rossler x 分量互信息曲线);画出来的曲线通常会在 tau8 到 tau12 之间出现第一个明显的低谷。保存的 tau 就是这个低谷对应的采样点数。注意这个 tau 是“采样点间隔数”不是时间单位。如果数据的 dt0.01那么实际延迟时间是 tau*dt0.08 左右。如果曲线没有清晰低谷而是整体平滑下掉最常见的原因是数据量太大导致分箱太粗或者序列本身不是混沌信号。把 bins 改成 24 再试一次同时缩短 maxTau 到 100往往就能看到局部极小。4.3 运行 cao_Single 扫描嵌入维数拿到 tau 之后把它写进脚本或直接在命令行调用maxDim 15; [E1, E2] cao_Single(data, tau, maxDim); figure; plot(1:maxDim, E1, o-, 1:maxDim, E2, s-); legend(E1, E2); xlabel(d);观察 E2 曲线从 d1 开始E2 一般大于 1随着 d 增加逐渐下降当它首次落到 1 附近并保持时这个 d 就是推荐的嵌入维数。对 Rossler x 分量通常出现在 3 到 5 之间。如果 E2 一直在 1.5 以上不掉说明 tau 取得太大回到 4.2 把 tau 减小一两个点重跑。4.4 解读 m1m2 的值.txt 与 m1.fig资源包里预留了 m1m2 的值.txt 和 m1.fig这是作者跑完 doubleCao.m 后留下的标准输出。用任何文本编辑器打开 m1m2 的值.txt典型内容是两个数字比如4 4第一行是 m1第二行是 m2。含义如下文件字段对应曲线判据注意事项m1E1 平台起点abs(diff(E1)) 0.05对数据噪声较敏感m2E2 稳定点abs(E2-1) 0.2推荐优先使用m1.fig 是 MATLAB 图窗文件用 openfig(m1.fig) 打开。如果打开后 E1、E2 曲线在 d 比较小时交叉得很乱那大概率是生成数据时瞬态没去掉回到 data_Rossler.m 把 x 的截取起点从 5000 调到 8000 再生成一次。5. 避坑嵌入维数计算中五个常见翻车现场5.1 采样步长 dt 没写进延迟把采样点数当时间现象算出 tau10直接当作延迟时间 10 秒送进重构函数结果相图拓扑完全不对Cao 方法给出的 m 也毫无规律。原因互信息计算输出的 tau 是采样点索引不是物理时间。Rossler 系统的习惯采样步长 dt0.01 时tau10 实际只对应 0.1 个时间单位。解决统一单位。延迟要么全程用采样点数要么全程用时间单位。我的习惯是在 reconstitution 内部只用采样点数所有对外沟通都换算成实际时间实际延迟 tau * dt。5.2 最近邻找成了自己E1 恒等于 1现象E1 曲线从头到尾都是 1几乎没有任何上升E2 也全部等于 1程序“完美”地给出错误结论。原因knnsearch 返回 K1 的最近邻而这个邻居就是状态点自身距离比为 1。这是 Cao 方法最常见的翻车点比阈值问题隐蔽得多。解决knnssearch(Xd, Xd, K, 2)取第 2 列作为邻居编号。还有一种写法是 idx knnsearch(Xd, Xd, K, 2); idx(i, 2) 里就是我们要的 j。用完后可以在命令行验证min(idx(:, 2)) 应该大于 0且大多数不等于 1。5.3 数据量太少E2 曲线剧烈抖动现象N1000 时E2 在 d3 附近掉到 1d7 又跳回 1.8整体像噪声。原因Cao 方法本质是统计均值N 太小导致状态点数量 num_points 太少最近邻距离比值的方差压不下去。数据量少于 5000 点时经常出现假收敛。解决把 data_Rossler.m 的 T 加长比如 T0:dt:500截取瞬态后保留 40000 点左右。如果处理的是实测信号想办法补长记录或者先用三次样条插值后再计算插值不会增加真实信息量但能让 E2 曲线稳定一些。5.4 reconstitution 的行列顺序写反程序直接报错现象运行 cao_Single 时报 Index exceeds the number of array elements断点定位在 X(i, k) data(i (k-1)*tau) 这一行。原因调用 reconstitution(data, d, tau) 时把 m 和 tau 的位置写反了导致 num_points 计算出来的值超过 data 长度。解决固定函数签名是 reconstitution(data, m, tau)不要学某些开源代码把 tau 放在 m 前面。如果已经反了最简单的修复是把调用处参数换回来或在 reconstitution 里加一个防御if nargin 3 tau m; m 2; % 或从参数名推断 end更推荐的做法是统一所有调用点不要依赖这种隐式修复。5.5 用自相关函数过零点当延迟Rossler 吸引子被压扁现象tau1 或 2重构出来的三维相图是一团靠近对角线的点看不出 Rossler 吸引子特有的带状螺旋。原因自相关函数只反映线性相关对 Rossler 这类非线性系统第一个过零点往往太小导致相邻坐标之间冗余度过高。解决改用互信息曲线第一个极小值。如果手头没有互信息函数就用 3.4 的版本。这个坑我在最初做混沌分析时踩过后来养成了“先互信息后 Cao”的固定流程两种方法结果差异最大时可以差到 2 个维数。6. 多一步验证用替代数据法确认嵌入维数不是巧合6.1 洗掉时间顺序看 E2 曲线是否变样Cao 方法给出来的 m 是一个建议值它有没有意义取决于数据里是否真的存在动力学结构。一个简单有效的验证办法是替代数据法把时间序列随机打乱破坏所有时间相关性然后对打乱后的序列重新跑一遍 Cao。rng(1); surrogate data(randperm(length(data))); [E1s, E2s] cao_Single(surrogate, tau, maxDim); figure; plot(1:maxDim, E2, o-, 1:maxDim, E2s, s-); legend(原始数据 E2, 替代数据 E2);正常情况是原始数据的 E2 会在某个 d 附近快速掉到 1 并保持替代数据的 E2 没有清晰转折要么一直波动要么下降得缓慢而不定型。如果两者长得很像说明你手里这段序列可能不是确定性混沌信号Cao 算出的维数没有物理意义。这一步不必追求统计检验目测曲线差异就能过滤掉一半以上的误判。6.2 用相图目测做最终确认除了替代数据我每次拿到 m 之后还会做一次“目检”。用 reconstitution 把序列重构到 Cao 建议的维数取前三个坐标画三维相图m 4; % 用 Cao 得到的值 X reconstitution(data, m, tau); plot3(X(:, 1), X(:, 2), X(:, 3), .); xlabel(x(t)); ylabel(x(ttau)); zlabel(x(t2tau));Rossler 系统的吸引子应该呈现清晰的单瓣螺旋边界锐利没有明显毛刺。如果相图像一坨云说明 m 或 tau 至少有一个不合适如果相图折叠严重说明 m 偏小。这个检查不需要任何额外工具箱肉眼几十秒就能判断。6.3 把 m 写进后续分析确认输出稳定最后一步很实际把 Cao 给出的 m 代入之后的 Lyapunov 指数计算或预测模型看结果有没有落在合理区间。Rossler 系统的最大 Lyapunov 指数通常在 0.07 到 0.13 之间如果算出来是负数或者大得离谱优先回头检查重构参数而不是怀疑算法。从那以后我每次跑 Cao 方法都强制走一遍“替代数据对照 相图目检 Lyapunov 合理性检查”这三步半个小时能避免整周的返工。希望帮到你。本文还有配套的精品资源点击获取 简介面向混沌时间序列分析研究者与MATLAB使用者提供基于Cao方法求取嵌入维数与时间延迟的完整实现聚焦Rossler系统的相空间重构与互信息计算。压缩包共9个文件以6个m脚本为核心覆盖数据生成、互信息计算、Cao算法主流程和相空间重构两套计算流程便于对照验证另附2个txt中间结果和1个fig图形文件用于核对输出整包仅8KB运行轻量、易于二次修改。已有356人学习下载。Cao方法通过互信息曲线极小值确定最佳延迟并结合嵌入维数判断准则完成状态空间重构脚本保留从Rossler时间序列生成、延迟量选取到维数验证的完整链路结果图可帮助理解参数变化对重构质量的影响。替换输入序列后可迁移至其他混沌动力学系统的延迟与维数估计适合科研入门、课程实验或复现经典混沌分析案例。1. 嵌入维数这道坎用 Cao 方法解 Rossler 系统的重构参数从一段混沌时间序列里恢复动力学信息第一步就卡在相空间重构上嵌入维数 m 和时间延迟 τ 取多少直接决定后面算 Lyapunov 指数、关联维数、做预测的结果。我用自相关法取过延迟也被 FNN 方法的阈值折腾过——阈值一改维数就变基本靠玄学。后来换成 Cao 方法配合互信息找延迟才把这组参数稳定住。这份 cao-method MATLAB 资源就是解决这个问题的它自带 Rossler 系统生成、相空间重构、Cao 维数扫描脚本输出 m1/m2 值适合做混沌动力学、非线性时间序列和信号重构方向的从业者。你不需要自己从头推导改改参数就能复现。2. 原理先立住互信息定延迟、Cao 法定维数为什么能避开阈值2.1 时间延迟不是越小越好互信息找第一个极小值相空间重构是把一维标量序列 x(t) 变成一组 m 维状态向量 [x(t), x(tτ), ..., x(t(m-1)τ)]目的是还原原系统的吸引子拓扑。这里 τ 的选择非常微妙τ 太小相邻坐标几乎一样吸引子被压扁在对角线附近重构出来的点都在一条线上τ 太大相邻坐标之间已经失去关联变成随机散射动力学信息被噪声盖住。自相关函数 R(τ) 是很多人第一个会试的工具通常取它第一次过零点的位置当延迟。问题在于自相关只刻画线性相关性而 Rossler 系统的动力学本质是非线性的线性无关不等于统计独立。我实际跑下来的体验是自相关给出的 τ 往往偏小重构出的 Rossler 吸引子明显挤压横幅比正常要扁。互信息是更适合非线性系统的工具。它的定义是I(τ) Σ p(x(t), x(tτ)) · log[p(x(t), x(tτ)) / (p(x(t)) · p(x(tτ)))]它衡量的是两个变量之间的统计依赖程度不限于线性关系。p(x(t)) 和 p(x(tτ)) 是一维概率p(x(t), x(tτ)) 是联合概率。实际计算时先把数据范围分成若干等宽箱子典型取 16 或 32统计数据点落在每个箱的频率再用这些频率近似概率。取 I(τ) 第一次出现局部极小值的 τ作为重构延迟。这个位置意味着在保留动态依赖的前提下坐标间的冗余信息最少。代码层面我一般把最大搜索范围设为总长度的 5%10%比如 20000 点就搜到 200 左右。箱子数不是越多越好数据量只有几千点时设 32 个箱子会让很多格子出现零频率估计出的互信息会抖动得很厉害。2.2 Cao 方法用 E1/E2 比值替代 FNN 阈值确定延迟 τ 之后剩下的嵌入维数 m 是第二个坎。最常用的 FNN假近邻法思路是在 m 维空间找每个点的最近邻然后看这两个点在 m1 维空间里是否还是近邻如果不再是说明原来那个邻居是因为维数不够而造成的“假邻居”。问题在于判断“还是不是近邻”需要给定一个相对阈值 R_tol 和绝对阈值 A_tol这两个阈值完全是主观的。阈值设 10 或设 15同一个序列能算出两个不同的维数。Cao 方法Cao, 1997就是为了消掉这个主观阈值。它的做法是这样对每个候选维数 d先把序列重构到 d 维和 d1 维然后在 d 维空间里找每个状态点的最近邻要排除自己接着计算这两个邻近点在 d1 维空间里的距离与它们在 d 维空间里距离的比值a(i,d) || X_i(d1) - X_j(d1) || / || X_i(d) - X_j(d) ||其中 j 是 i 的最近邻索引。对所有状态点取平均得到 E1(d)E1(d) (1/(N - dτ)) · Σ a(i,d)当 d 小于真实嵌入维数时维度每增加一维原来被压扁的近邻会被拉开E1 会明显变化当 d 已经达到真实维数之后再增加维度最近邻的相对距离不再系统性变大E1 进入平台。更常用的判据是 E2E2(d) E1(d1) / E1(d)当 E2 从明显大于 1 的状态掉到 1 附近并保持稳定就说明已经到达真实嵌入维数附近。Cao 方法不需要设阈值只看曲线趋势这就是它在工程里比 FNN 顺手的原因。需要注意Cao 方法输出的是“最大需要的维数”而不是“系统真实维数”。Rossler 系统的吸引子分形维数大概在 2.1 左右但 Takens 嵌入定理要求 m ≥ 2D1所以理论下限是 5 到 6实际重构时因为数据有限Cao 给出的 m 通常在 36 之间。你最后用的是 m 而不是 D这一点别混。2.3 为什么选 Rossler 当测试对象以及 m 与 τ 的连带关系Rossler 系统是经典三维混沌模型方程是 dx/dt -y-zdy/dt xa·ydz/dt bz·(x-c)取 a0.2、b0.2、c5.7 时处在混沌状态。它比 Lorenz 系统简单吸引子形态清楚也更容易目测验证所以最适合当 Cao 方法的测试对象。这份资源里所有脚本都围绕这条主线设计你可以先用它跑通流程、看懂输出再把自己的数据替换到 data 变量里。还有一点容易被忽略延迟 τ 和嵌入维数 m 不是两个独立变量。τ 取小的时候Cao 的 E2 曲线往往更早进入平台给人“维数很小”的错觉τ 取太大数据点之间相关性太弱E1 曲线上会出现不必要的波动。所以我的习惯是先算互信息把 τ 固定下来再跑 Cao如果算出 m 之后发现相图挤压或发散再回头微调 τ 扫一次看 m 是否变化。这不是严谨的理论结论但实际处理混沌数据时很实用。3. 拆 MATLAB 源码data_Rossler、reconstitution、cao_Single 在干什么3.1 data_Rossler.m先解决“数据从哪来”资源包里的 data_Rossler.m 是入口脚本负责生成用于验证的 Rossler 时间序列。直接用 MATLAB 的 ode45 积分即可关键是把输出处理成干净的标量时间序列。% data_Rossler.m % 生成 Rossler 系统的 x 分量时间序列 clear; clc; a 0.2; b 0.2; c 5.7; % Rossler 混沌参数 dt 0.01; % 采样步长 T 0:dt:200; % 总仿真时长 200 个时间单位 x0 [0; 0; 1]; % 初始状态 % 用 ode45 积分 Rossler 方程 [t, x] ode45((t,x) [ -x(2)-x(3); x(1)a*x(2); bx(3)*(x(1)-c) ], T, x0); % 去掉前 5000 点瞬态取 x 分量作为观测序列 data x(5000:end, 1); save(rossler_data.mat, data, dt);这段代码的逻辑不复杂先定义参数再用 ode45 积分最后保存数据。需要注意两个细节。第一c5.7 是 Rossler 系统进入混沌的经典取值不少人把 c 改成 4 发现系统进入周期窗口然后误以为 Cao 方法算错了维数。第二前 5000 点必须丢弃因为初值 [0;0;1] 到吸引子上的轨迹有一段瞬态过程保留它会让互信息和 E1 曲线都带上“污染段”。这里的 dt 不是数值积分的内部步长ode45 是自适应步长dt 只决定输出点的间隔。对 Rossler 系统dt0.01 已经足够密每个振荡周期大约能采 60100 个点。如果你把 dt 改到 0.05互信息曲线会明显变粗糙有可能错过第一个极小值点。采样总时长 200 个时间单位对应 20000 个数据点去瞬态后剩 15001 点这个量对 Cao 方法来说是合格的。3.2 reconstitution.m相空间重构核心函数reconstitution.m 是整个包的地基Cao 方法的每一步都要调用它。它的输入是原始序列、嵌入维数 m、延迟 tau输出是重构后的相空间矩阵。function X reconstitution(data, m, tau) % data: 原始标量序列列向量 % m: 嵌入维数 % tau: 延迟时间 % X: 重构后的相空间矩阵每一行是一个状态点 N length(data); num_points N - (m - 1) * tau; % 能构造出的状态点数量 X zeros(num_points, m); for i 1:num_points for k 1:m X(i, k) data(i (k - 1) * tau); end end end逻辑很简单num_points 是最多能取到的状态点数因为最后一个状态点要取到 data(N)所以索引是 i(m-1)τ ≤ N。循环里 i 表示状态点序号k 表示当前点在这一行的第几列。输出矩阵的行数是 num_points列数是 m这也是后续 knnsearch 能直接使用的格式。我见过不少人在这个函数上报错典型原因是把 num_points 和 m 写反。这里再强调一次XO reconstitution(data, m, tau) 之后size(X) 第一维是 N-(m-1)τ 而不是 m。后面做 Cao 计算时如果发现索引越界或者 mean(a) 结果异常请先回头打印一下 size(X)。3.3 cao_Single.m 与 doubleCao.m从单次扫描到全维曲线cao_Single.m 是核心计算函数它对每个候选维数 d 计算 E1 和 E2。我按资源包的名字和输出文件推断作者把“doubleCao”理解成同时计算 E1、E2 两条曲线的完整流程。这里给出一个能直接跑的兼容版本function [E1, E2] cao_Single(data, tau, maxDim) % 对每个嵌入维数 d 计算 E1(d) 和 E2(d) N length(data); E1 zeros(1, maxDim); for d 1:maxDim Xd reconstitution(data, d, tau); % d 维相空间 Xd1 reconstitution(data, d 1, tau); % d1 维相空间 % 找每个点的最近邻K2 时第 1 列是自己第 2 列才是真邻居 idx knnsearch(Xd, Xd, K, 2); sum_a 0; numP size(Xd, 1); for i 1:numP j idx(i, 2); dist_d norm(Xd(i, :) - Xd(j, :)); dist_d1 norm(Xd1(i, :) - Xd1(j, :)); sum_a sum_a dist_d1 / (dist_d eps); end E1(d) sum_a / numP; end E2(1) 1; for d 2:maxDim E2(d) E1(d) / E1(d - 1); end end逻辑说明对每个 d先在 d 维空间找最近邻 j再把这个邻居索引套到 d1 维空间里算距离比。这里的关键是 K2因为 K1 时最近邻就是自己距离比为 1E1 曲线会变成一根平线看起来像已经收敛实际是假象。注意 E2 的定义是 E1(d1)/E1(d)我在代码里用后一个 E1 除以前一个 E1 来填充 E2(d)方向是一致的。判断规则是当 E2 从明显大于 1 的位置掉下来然后在 1 附近持续波动这个转折点对应的 d 就是 Cao 方法建议的嵌入维数。doubleCao.m 是调用 cao_Single 并输出结果的上层脚本。原包里应该包含类似下面的逻辑只是文件命名和保存格式可能有差异% doubleCao.m 求 Rossler 系统嵌入维数 clear; clc; load(rossler_data.mat); tau 8; % 由互信息曲线第一个极小值得到 maxDim 15; [E1, E2] cao_Single(data, tau, maxDim); % 找 E1 进入平台的起始维数常见判据变化幅度小于 0.05 m1 find(abs(diff(E1)) 0.05, 1); % 找 E2 首次稳定在 1 附近的维数允许 0.2 的容差 m2 find(abs(E2 - 1) 0.2, 1); fprintf(m1 %d, m2 %d\n, m1, m2); % 保存结果 save(m1m2的值.txt, m1, m2, -ascii); save(m1的值.txt, m1, -ascii); % 画图并保存 figure; plot(1:maxDim, E1, o-, 1:maxDim, E2, s-); legend(E1, E2); xlabel(嵌入维数 d); ylabel(E1 / E2); savefig(m1.fig);m1 对应 E1 曲线进入平台的起点m2 对应 E2 稳定在 1 附近的起点。两者通常一致如果差超过 1优先信 m2因为 E1 平台的“判稳”比 E2 更依赖容差设置。你在资源包里看到的 m1.fig 就是这个图m1m2 的值.txt 是这两个数字。3.4 包里没有互信息函数时自己补一份文件列表里没有专门的互信息计算脚本但标题里写了“延迟 互信息”说明作者默认这一步已经单独做过了。我在这里补一个直方图分箱的简化实现方便你把流程串起来。function [tau, mi] mutual_info_delay(data, maxTau, bins) % 用等宽直方图分箱估计互信息返回曲线和第一个极小值 data data(:); N length(data); edges linspace(min(data), max(data), bins 1); dataIdx discretize(data, edges); numObs N - maxTau; mi zeros(1, maxTau); for tau 1:maxTau pxy zeros(bins, bins); for i 1:numObs pxy(dataIdx(i), dataIdx(i tau)) ... pxy(dataIdx(i), dataIdx(i tau)) 1; end pxy pxy / numObs; px sum(pxy, 2); % bins x 1 py sum(pxy, 1); % 1 x bins logTerm log(pxy ./ (px * py) eps); mi(tau) sum(pxy(:) .* logTerm(:)); end % 找第一个局部极小值数据质量差时可直接用全局最小 [~, tau] min(mi); end逻辑说明先对数据分箱再用二维频数近似联合概率 pxy一维频数近似边缘概率 px 和 py。px * py 是外积得到 bins×bins 的独立假设概率矩阵两者比值的对数就是互信息贡献。最后相加得到 I(τ)。参数上bins 通常取 16如果你的数据只有 3000 点用 8 到 10 个箱子更稳。maxTau 我一般取 N/20 左右太大会让 numObs 太小尾部统计不可靠。代码里取的是全局最小值严格做法应该找第一个局部极小点两者在 Rossler 这类清晰系统上差别不大。4. 实操全流程从生成 Rossler 序列到读出 m1/m2 值4.1 搭建目录并确认数据量拿到资源包后先把所有 .m 文件放进同一个目录MATLAB 里 cd 到该目录。然后运行 data_Rossler.m观察工作区里 data 的长度。这里有一个检查数据的习惯不要急着跑后续脚本先确认 data 是列向量、长度在 10000 以上。如果长度只有几百后面互信息和 Cao 都会给出离谱结果。cd(你的路径/cao-method); run(data_Rossler.m); whos data dt length(data) % 应该输出 15001运行后 rossler_data.mat 会出现在当前目录。这个文件就是后续所有计算的唯一数据源。如果你要处理自己的数据直接在这一步之后把 data 替换掉并把 dt 改成你自己的采样间隔。4.2 用互信息确定延迟 tau按照 3.4 的函数在命令行执行[tau, mi] mutual_info_delay(data, 200, 16); figure; plot(1:200, mi); xlabel(tau (采样点数)); ylabel(I(tau)); title(Rossler x 分量互信息曲线);画出来的曲线通常会在 tau8 到 tau12 之间出现第一个明显的低谷。保存的 tau 就是这个低谷对应的采样点数。注意这个 tau 是“采样点间隔数”不是时间单位。如果数据的 dt0.01那么实际延迟时间是 tau*dt0.08 左右。如果曲线没有清晰低谷而是整体平滑下掉最常见的原因是数据量太大导致分箱太粗或者序列本身不是混沌信号。把 bins 改成 24 再试一次同时缩短 maxTau 到 100往往就能看到局部极小。4.3 运行 cao_Single 扫描嵌入维数拿到 tau 之后把它写进脚本或直接在命令行调用maxDim 15; [E1, E2] cao_Single(data, tau, maxDim); figure; plot(1:maxDim, E1, o-, 1:maxDim, E2, s-); legend(E1, E2); xlabel(d);观察 E2 曲线从 d1 开始E2 一般大于 1随着 d 增加逐渐下降当它首次落到 1 附近并保持时这个 d 就是推荐的嵌入维数。对 Rossler x 分量通常出现在 3 到 5 之间。如果 E2 一直在 1.5 以上不掉说明 tau 取得太大回到 4.2 把 tau 减小一两个点重跑。4.4 解读 m1m2 的值.txt 与 m1.fig资源包里预留了 m1m2 的值.txt 和 m1.fig这是作者跑完 doubleCao.m 后留下的标准输出。用任何文本编辑器打开 m1m2 的值.txt典型内容是两个数字比如4 4第一行是 m1第二行是 m2。含义如下文件字段对应曲线判据注意事项m1E1 平台起点abs(diff(E1)) 0.05对数据噪声较敏感m2E2 稳定点abs(E2-1) 0.2推荐优先使用m1.fig 是 MATLAB 图窗文件用 openfig(m1.fig) 打开。如果打开后 E1、E2 曲线在 d 比较小时交叉得很乱那大概率是生成数据时瞬态没去掉回到 data_Rossler.m 把 x 的截取起点从 5000 调到 8000 再生成一次。5. 避坑嵌入维数计算中五个常见翻车现场5.1 采样步长 dt 没写进延迟把采样点数当时间现象算出 tau10直接当作延迟时间 10 秒送进重构函数结果相图拓扑完全不对Cao 方法给出的 m 也毫无规律。原因互信息计算输出的 tau 是采样点索引不是物理时间。Rossler 系统的习惯采样步长 dt0.01 时tau10 实际只对应 0.1 个时间单位。解决统一单位。延迟要么全程用采样点数要么全程用时间单位。我的习惯是在 reconstitution 内部只用采样点数所有对外沟通都换算成实际时间实际延迟 tau * dt。5.2 最近邻找成了自己E1 恒等于 1现象E1 曲线从头到尾都是 1几乎没有任何上升E2 也全部等于 1程序“完美”地给出错误结论。原因knnsearch 返回 K1 的最近邻而这个邻居就是状态点自身距离比为 1。这是 Cao 方法最常见的翻车点比阈值问题隐蔽得多。解决knnssearch(Xd, Xd, K, 2)取第 2 列作为邻居编号。还有一种写法是 idx knnsearch(Xd, Xd, K, 2); idx(i, 2) 里就是我们要的 j。用完后可以在命令行验证min(idx(:, 2)) 应该大于 0且大多数不等于 1。5.3 数据量太少E2 曲线剧烈抖动现象N1000 时E2 在 d3 附近掉到 1d7 又跳回 1.8整体像噪声。原因Cao 方法本质是统计均值N 太小导致状态点数量 num_points 太少最近邻距离比值的方差压不下去。数据量少于 5000 点时经常出现假收敛。解决把 data_Rossler.m 的 T 加长比如 T0:dt:500截取瞬态后保留 40000 点左右。如果处理的是实测信号想办法补长记录或者先用三次样条插值后再计算插值不会增加真实信息量但能让 E2 曲线稳定一些。5.4 reconstitution 的行列顺序写反程序直接报错现象运行 cao_Single 时报 Index exceeds the number of array elements断点定位在 X(i, k) data(i (k-1)*tau) 这一行。原因调用 reconstitution(data, d, tau) 时把 m 和 tau 的位置写反了导致 num_points 计算出来的值超过 data 长度。解决固定函数签名是 reconstitution(data, m, tau)不要学某些开源代码把 tau 放在 m 前面。如果已经反了最简单的修复是把调用处参数换回来或在 reconstitution 里加一个防御if nargin 3 tau m; m 2; % 或从参数名推断 end更推荐的做法是统一所有调用点不要依赖这种隐式修复。5.5 用自相关函数过零点当延迟Rossler 吸引子被压扁现象tau1 或 2重构出来的三维相图是一团靠近对角线的点看不出 Rossler 吸引子特有的带状螺旋。原因自相关函数只反映线性相关对 Rossler 这类非线性系统第一个过零点往往太小导致相邻坐标之间冗余度过高。解决改用互信息曲线第一个极小值。如果手头没有互信息函数就用 3.4 的版本。这个坑我在最初做混沌分析时踩过后来养成了“先互信息后 Cao”的固定流程两种方法结果差异最大时可以差到 2 个维数。6. 多一步验证用替代数据法确认嵌入维数不是巧合6.1 洗掉时间顺序看 E2 曲线是否变样Cao 方法给出来的 m 是一个建议值它有没有意义取决于数据里是否真的存在动力学结构。一个简单有效的验证办法是替代数据法把时间序列随机打乱破坏所有时间相关性然后对打乱后的序列重新跑一遍 Cao。rng(1); surrogate data(randperm(length(data))); [E1s, E2s] cao_Single(surrogate, tau, maxDim); figure; plot(1:maxDim, E2, o-, 1:maxDim, E2s, s-); legend(原始数据 E2, 替代数据 E2);正常情况是原始数据的 E2 会在某个 d 附近快速掉到 1 并保持替代数据的 E2 没有清晰转折要么一直波动要么下降得缓慢而不定型。如果两者长得很像说明你手里这段序列可能不是确定性混沌信号Cao 算出的维数没有物理意义。这一步不必追求统计检验目测曲线差异就能过滤掉一半以上的误判。6.2 用相图目测做最终确认除了替代数据我每次拿到 m 之后还会做一次“目检”。用 reconstitution 把序列重构到 Cao 建议的维数取前三个坐标画三维相图m 4; % 用 Cao 得到的值 X reconstitution(data, m, tau); plot3(X(:, 1), X(:, 2), X(:, 3), .); xlabel(x(t)); ylabel(x(ttau)); zlabel(x(t2tau));Rossler 系统的吸引子应该呈现清晰的单瓣螺旋边界锐利没有明显毛刺。如果相图像一坨云说明 m 或 tau 至少有一个不合适如果相图折叠严重说明 m 偏小。这个检查不需要任何额外工具箱肉眼几十秒就能判断。6.3 把 m 写进后续分析确认输出稳定最后一步很实际把 Cao 给出的 m 代入之后的 Lyapunov 指数计算或预测模型看结果有没有落在合理区间。Rossler 系统的最大 Lyapunov 指数通常在 0.07 到 0.13 之间如果算出来是负数或者大得离谱优先回头检查重构参数而不是怀疑算法。从那以后我每次跑 Cao 方法都强制走一遍“替代数据对照 相图目检 Lyapunov 合理性检查”这三步半个小时能避免整周的返工。希望帮到你。本文还有配套的精品资源点击获取