简介：面向混沌时间序列分析与非线性动力学研究，这份MATLAB资源基于Cao方法实现嵌入维数与延迟时间的自动估计，并以经典Rossler系统作为验证对象。资源包共9个文件，总大小仅8KB，其中6个m脚本覆盖数据生成、互信息计算、Cao算法主程序与参数求解，2个txt文件保存计算所得的m1/m1m2数值，1个fig图展示嵌入维数曲线，结构紧凑、便于对照运行。目前已有355人学习。通过调用这套代码，读者可直观理解Cao方法的完整流程，快速获得Rossler系统的最佳时间延迟与嵌入维数，同时还能替换数据生成模块以拓展至其他混沌系统，适合需要掌握相空间重构技术的学生和科研人员参考。
1. 嵌入维数估计为什么不能只靠眼睛——Cao 方法与互信息的定位
拿到一段 Rossler 系统的 x 分量时间序列，第一反应通常是直接画相图、看轨迹，然后凭经验猜一个嵌入维数。这种做法在噪声小、数据长的实验室数据上勉强能看，但一旦换到实测信号，比如生物电信号或者金融收益率序列，同一个系统在不同时间段里画出来的重构轨迹可能完全不像同一个东西。问题不在于系统变了，而在于延迟 τ 和嵌入维数 m 没选对。Cao 方法（Cao's method）解决的就是这个问题：它用一组连续变化的量 E1(d) 和 E2(d) 来判断嵌入维数什么时候饱和，不需要人为设定阈值；而互信息法（mutual information）则负责在重构之前先把延迟 τ 定下来，它的核心是找互信息曲线上的第一个极小值。这套组合拳在混沌时间序列分析里几乎是标准预处理流程，对做混沌控制、预测模型输入维度选择、以及动力学不变量估计的人都适用。下面从数据生成开始，逐步拆解这个压缩包里每个脚本的实际作用。
2. Rossler 系统数据生成与相空间重构：data_Rossler.m 和 reconstitution.m 在做什么
Cao 方法的输入不是原始微分方程，而是采样后的一维时间序列。所以第一步必须先把 Rossler 系统的连续微分方程离散化，得到足够长的观测数据，再通过延迟坐标重构把一维序列映射到高维空间。这一步做不好，后面算出来的互信息和嵌入维数都会失真。
2.1 从 ODE 到时间序列：data_Rossler.m 的采样与参数设置
Rossler 系统是最经典的混沌系统之一，它的状态方程只有三个变量，但动力学行为比 Lorentz 系统更容易调节。标准形式是：
dx/dt = -y - z dy/dt = x + a*y dz/dt = b + z*(x - c)经典参数取 a = 0.2，b = 0.2，c = 5.7，此时系统处于混沌状态。data_Rossler.m 这个脚本的核心任务就是利用 MATLAB 的 ode45 求解这个方程组，然后按固定采样间隔抽取 x 分量作为后续分析的时间序列。
% data_Rossler.m % 生成Rossler系统x分量时间序列，用于后续Cao方法与互信息计算 clear; clc; % 系统参数：a=0.2, b=0.2, c=5.7 对应混沌状态 a = 0.2; b = 0.2; c = 5.7; % 定义Rossler微分方程 rossler = @(t, x) [-x(2) - x(3); x(1) + a * x(2); b + x(3) * (x(1) - c)]; % 初始条件与时间跨度 x0 = [0; 0; 0]; % 初始状态 t_total = 200; % 总仿真时长 [t, y] = ode45(rossler, [0 t_total], x0); % 数值积分 % 采样：丢掉前50个时间单位，消除暂态 skip = 50; idx = find(t > skip); data = y(idx, 1); % 取x分量 % 降采样：等间隔抽取，避免相邻点过密导致互信息计算失真 fs = 10; % 采样频率 data = data(1:fs:end); save('Rossler_x.mat', 'data');代码的逻辑很直接：先用ode45做数值积分得到连续轨迹，然后通过find(t > skip)去掉初始暂态部分。跳过暂态这一步很关键，因为从 (0,0,0) 出发的轨迹还没有落到吸引子上，前面这一段不属于稳态混沌运动，直接保留会污染延迟和嵌入维数的估计。
降采样的参数fs = 10表示每 10 个积分步取一个点。这里需要说明一个常见误区：ode45的步长是自适应变化的，如果直接把y全部拿出来用，相邻两个点的时间间隔不均匀，互信息法假设等间隔采样，不均匀序列会让延迟 τ 的物理意义变得模糊。所以我在降采样时用的是每隔固定点数抽取，实际使用时如果原始积分步长变化大，更稳妥的做法是先用interp1重采样到固定时间网格，再做抽取。
2.2 reconstitution.m：延迟坐标重构的矩阵构造
有了时间序列data，下一步就是构造延迟坐标向量。Takens 嵌入定理说，如果原始动力学系统的吸引子维数是 d，那么嵌入维数 m 只要满足 m > 2d 就能把吸引子拓扑还原出来。重构方式是把一维序列x(1), x(2), ..., x(N)变成一组 m 维向量：
X(i) = [x(i), x(i + τ), x(i + 2τ), ..., x(i + (m-1)τ)]reconstitution.m 就是做这件事的。它的输入是原始序列、延迟 τ 和嵌入维数 m，输出是重构后的相空间矩阵，每一行是一个相点。
% reconstitution.m % 延迟坐标重构：把一维时间序列映射为m维相空间 function X = reconstitution(data, tau, m) N = length(data); % 有效相点数量：末尾不足(m-1)*tau的部分丢弃 M = N - (m - 1) * tau; X = zeros(M, m); for i = 1:M for j = 1:m X(i, j) = data(i + (j - 1) * tau); end end end这个函数实现的是标准的 Takens 重构。参数tau是延迟步数，m是嵌入维数，M是重构后相点的个数。需要注意M = N - (m - 1) * tau这个关系：序列开头和末尾都有无法完整构造相点的区域，开头丢掉的是 0 到 τ 之间的索引，末尾丢掉的是最后(m-1)*tau个样本。这也是为什么嵌入维数不能取得过大——m 越大，有效相点越少，后续距离计算和邻居搜索的统计基础就越弱。
实际使用时，重构矩阵 X 的每一行就代表高维空间中的一个状态点。Cao 方法后续计算欧氏距离、判断最近邻居，都是在 X 这个矩阵上做的，它本质上决定了后续所有统计量的可靠性。
3. 用互信息曲线找延迟 τ：第一极小值判据与 MATLAB 实现
延迟 τ 的选择直接影响重构质量。如果 τ 太小，相邻延迟坐标之间几乎完全相关，重构吸引子被压缩在对角线附近；如果 τ 太大，混沌系统的相邻状态在延迟 τ 后已经指数分离，重构出的相点之间看起来像是随机噪声。互信息法通过衡量原始序列与延迟序列之间的统计依赖程度来找到这两者之间的平衡点。
3.1 互信息的定义与直方图估计
互信息是信息论里量化两个随机变量相关性的指标，和线性相关性不同，它不假设变量之间的关系是线性的。这一点在混沌时间序列里特别重要，因为混沌吸引子上的变量关系天然是非线性的，用线性自相关函数算出来的延迟往往偏大或偏小。
互信息的定义是：
I(τ) = Σ P(x(i), x(i+τ)) * log( P(x(i), x(i+τ)) / (P(x(i)) * P(x(i+τ))) )其中联合概率P(x(i), x(i+τ))表示原始序列和延迟序列同时取某对值的概率。实际计算中概率分布未知，常见的做法是把信号分成多个区间（bin），统计每个区间内的点数，用频率近似概率。MATLAB 里可以用histogram或histcounts2做二维直方图：
% mutual_info.m 核心计算逻辑 % 输入：时间序列data，最大延迟max_tau，分箱数nbins % 输出：不同tau下的互信息I function I = mutual_info(data, max_tau, nbins) data = data(:); N = length(data); I = zeros(1, max_tau); % 归一化到[0,1]区间，方便统一分箱 dmin = min(data); dmax = max(data); data_norm = (data - dmin) / (dmax - dmin) * nbins + 0.5; data_norm = floor(data_norm); % 每个点落到某个bin for tau = 1:max_tau n = N - tau; x = data_norm(1:n); y = data_norm(1 + tau:n + tau); % 二维直方图：联合概率 joint = histcounts2(x, y, 0:nbins, 0:nbins); joint = joint / sum(joint(:)); % 归一化为联合分布 % 计算边缘概率 px = sum(joint, 2); py = sum(joint, 1); % 互信息：对每个非零概率项累加 [ix, iy] = meshgrid(1:nbins, 1:nbins); valid = joint > 0; I(tau) = sum(joint(valid) .* log(joint(valid) ./ (px(ix(valid)) .* py(iy(valid))))); end end代码里最关键的是histcounts2这一步，它同时统计了(x(i), x(i+τ))落在二维网格每个格子里的次数。nbins的选择直接影响结果：bin 数太少，概率估计失真，曲线过于平滑；bin 数太多，每个格子里样本点太少，统计涨落变大，曲线出现大量假极小值。对于 Rossler 系统这样的低维混沌系统，数据量在一万点左右时，nbins = 16或32是比较合理的起点。
3.2 不同 τ 下的互信息曲线与第一极小值
计算好I(τ)后，绘制曲线并找到第一个极小值点。这里的逻辑是：当 τ 很小时，x(i)和x(i+τ)几乎相同，互信息很大；随着 τ 增加，两者间相关性下降，互信息减小；到达第一个极小值后，由于混沌系统的轨道折叠特性，互信息可能出现局部反弹。取第一个极小值作为延迟，是为了在信息保留最大和去冗余之间找到平衡——如果跳过第一个极小值去取全局最小，τ 可能过大，导致相邻延迟坐标之间的实际关联已经丢失。
% 调用示例：计算互信息并定位第一极小值 max_tau = 100; nbins = 32; I = mutual_info(data, max_tau, nbins); % 绘图 figure; plot(1:max_tau, I, 'b-', 'LineWidth', 1.2); xlabel('\tau'); ylabel('互信息 I(\tau)'); title('Rossler x 分量互信息曲线'); % 找第一极小值：排除tau=1后，找第一个局部极小 I_smooth = smoothdata(I, 'gaussian', 5); % 平滑去毛刺 dI = diff(I_smooth); first_min_idx = find(dI(1:end-1) < 0 & dI(2:end) > 0, 1, 'first') + 1; fprintf('最佳延迟 tau = %d\n', first_min_idx);平滑是必要的。直接对原始互信息曲线求一阶差分找极小值，容易因为统计涨落定位到 τ = 2 或 τ = 3 这样的假极小点。我用smoothdata做一个高斯窗口平滑，窗口宽度 5 对曲线形状影响不大，但能滤掉大部分高频抖动。
3.3 延迟 τ 对后续 Cao 方法的影响
有个细节容易被忽略：Cao 方法里的 E1(d) 和 E2(d) 计算重度依赖最近邻距离，而最近邻是在延迟坐标构成的相空间里搜索的。如果 τ 选得太小，重构相空间中的点在取范数时各维度的值高度接近，最近邻搜索几乎等效于在一维直线上找邻居，有效维度不足；如果 τ 选得过大，由于混沌吸引子的有界性，延迟坐标之间出现大量伪交点，最近邻距离会被系统性拉大，导致 E1(d) 的收敛速度变慢。所以互信息和 Cao 方法不是各算各的，而是前后衔接的两段式流程。
提示：实际工程里如果后续要做 Lyapunov 指数估计，τ 的宽容度比做预测模型要大。预测模型对 τ 敏感，因为输入特征的冗余度直接影响回归问题的条件数。
4. Cao 方法求嵌入维数：cao_Single.m 与 doubleCao.m 的计算流程
Cao 方法的核心思想是看「增加嵌入维度时，相空间中的距离结构是否发生根本变化」。它不需要像虚假最近邻点法（FNN）那样设定距离阈值，所以对噪声的适应能力更强，在实测信号上更稳健。
4.1 为什么不用虚假最近邻点法
FNN 的原理是：如果嵌入维数不够，两个在高维空间中相距很远的点会因为低维投影而成为最近的邻居，增加维数后这些「虚假邻居」会突然分开。判断假邻居需要设定一个距离比阈值，通常取 10 或 15，但这个阈值在噪声环境下很不稳定——信噪比变化时同样一组数据的 FNN 比例曲线会整体平移，阈值固定后得到的嵌入维数要么偏大要么偏小。Cao 方法规避了这个主观因素，它定义了两个无量纲量 E1(d) 和 E2(d)：
E1(d) = a(i, d) 的均值，其中 a(i, d) = ||Y_i(d+1) - Y_neighbor(i)(d+1)|| / ||Y_i(d) - Y_neighbor(i)(d)||当 d 达到某个值后，E1(d) 趋于平稳，这个值就是嵌入维数。E2(d) 则用来区分确定性混沌信号和随机信号，如果 E2(d) 在某个 d 之后不趋于 1，说明信号是确定性的。
4.2 cao_Single.m 的核心实现
压缩包里的 cao_Single.m 负责在给定延迟 τ 的情况下计算 E1 和 E2 随嵌入维数 d 的变化曲线。
% cao_Single.m % 计算单个延迟tau下的E1(d)与E2(d)，用于确定嵌入维数 function [E1, E2] = cao_Single(data, tau, max_d) N = length(data); E1 = zeros(1, max_d); E2 = zeros(1, max_d); for d = 1:max_d % 构造 d 维和 d+1 维的延迟向量 Xd = reconstitution(data, tau, d); Xd1 = reconstitution(data, tau, d + 1); M = size(Xd, 1); % 对每个向量找最近邻（在d维空间） a_sum = 0; e_sum = 0; for i = 1:M % 计算Xd第i行到其余所有行的欧氏距离 diff = Xd - Xd(i, :); dist = sqrt(sum(diff.^2, 2)); % 排除自身：把自身距离设为无穷大 dist(i) = inf; % 找最近邻索引 [~, n_idx] = min(dist); % 计算a(i,d)，分子用d+1维空间的距离 dist_d = dist(n_idx); diff_d1 = Xd1(i, :) - Xd1(n_idx, :); dist_d1 = sqrt(sum(diff_d1.^2, 2)); if dist_d > 0 a_i = dist_d1 / dist_d; else a_i = inf; % 完全重合点，跳过或标记 end a_sum = a_sum + a_i; % E2用相邻距离a_i的比值累积，公式见Cao论文 e_sum = e_sum + abs(a_i - mean_a); % 累计偏差项（示例） end E1(d) = a_sum / M; E2(d) = e_sum / M; end end这里有几个实现时必须注意的点。第一，dist(i) = inf排除自身是必须的，否则每个点的最近邻都是自己，距离恒为 0，E1 失去意义。第二，分子分母都用欧氏距离，但因为 Xd1 比 Xd 多一列，所以分子天然比分母大，E1 的值普遍大于 1，这不影响判断饱和趋势。第三，mean_a需要维护一个滑动均值，实际代码里更常用的版本是先算出所有a_i存储在数组里，再统一计算均值和方差，上面这段是示范循环结构。
真正检验嵌入维数是否收敛的标准是观察 E1 曲线：从 d = 1 开始，E1 快速下降，当 d 达到某个值后，E1 的波动进入一个微弱的缓变带。这个缓变带开始的 d 就是嵌入维数。
4.3 doubleCao.m 与 m1m2.txt 的生成逻辑
包里还有个 doubleCao.m，这个脚本把整个流程串起来了。从文件名的 m1m2 可以看出，它要对多组延迟 τ 分别跑 Cao 方法，然后把每组得到的最佳嵌入维数 m1、m2 输出成文本文件。
% doubleCao.m % 对比不同tau下的Cao法结果，输出m1和m2 clear; clc; load('Rossler_x.mat', 'data'); tau_range = [5, 8, 10, 12, 15]; % 候选延迟列表 max_d = 12; results = []; for k = 1:length(tau_range) tau = tau_range(k); [E1, E2] = cao_Single(data, tau, max_d); % 找E1饱和点：从第2个点开始，连续3次变化小于2%判定饱和 m_opt = 1; for d = 3:max_d chg = abs(E1(d) - E1(d-1)) / E1(d-1); if chg < 0.02 m_opt = d - 1; break; end end results = [results; tau, m_opt]; end % 输出到文本 fid = fopen('m1m2的值.txt', 'w'); fprintf(fid, 'tau\tm1\tm2\n'); for k = 1:size(results, 1) fprintf(fid, '%d\t%d\t%d\n', results(k, 1), results(k, 2), results(k, 2)); end fclose(fid); disp('Cao方法计算结果已写入 m1m2的值.txt');判定饱和的这个 2% 阈值是经验值。Cao 原始论文里建议用 E1 曲线的形状变化来判断，没有给明确的数值标准。我一般会同时绘制 E1 曲线，用人工确认自动判定的结果，尤其是在数据长度较短或者噪声偏大的场景下，2% 的阈值可能提前判饱和。m1m2的值.txt里同时记录 m1 和 m2，是因为有时候第一个饱和点出现得比较早，需要用 E2 曲线的行为做二次确认——如果 E2 在 m1 处仍然不收敛到 1 附近，说明存在长程相关性，需要检查是否是数据段长度不够或者 τ 偏小。
| 参数 | 含义 | 推荐范围 | 调试建议 |
|---|---|---|---|
| tau | 延迟步数 | 由互信息第一极小值确定 | τ 过小时 Cao 方法收敛慢，过大时 E1 曲线抖动加剧 |
| max_d | 最大嵌入维数 | 数据维度的 3～5 倍 | Rossler 系统取 10～15 足够 |
| 饱和判据 | E1 相对变化率 | 1%～3% | 噪声大时提高阈值，避免过度延迟 |
| 数据长度 N | 序列长度 | 大于 5000 点 | 过短时最近邻统计不稳定，E1 曲线毛刺多 |
5. 两个高价值技巧：E1 曲线判读与典型错误规避
最后这部分不打算复述流程，而是分享两个实际调试时最值得留意的判断技巧。
5.1 看 E1 饱和位置时，顺带记录 E2 的值
很多人在用 Cao 方法的时候只记录 E1 开始饱和的 d，然后直接拿去重构。但 E2 曲线其实包含了重要的区分信息——E2(d) 在确定性混沌系统里不会全都等于 1，而是在某些 d 处显著偏离。如果 E2 在没有任何饱和迹象的区间里突然跳到接近 1 的位置，这通常说明数据里混入了强噪声。此时即使 E1 看起来收敛了，重构出的相空间大概率是噪声主导的。所以拿到 cao_Single.m 的输出后，我习惯把 E1 和 E2 画在同一张图的两个 subplot 里，以 E2 的偏离行为作为 E1 判据的交叉验证，而不是只看单一曲线。
5.2 三个容易踩的坑
第一个坑是直接用 ode45 的原始输出做互信息计算。自适应步长会让时间序列的相邻点间隔不相等，互信息曲线会出现周期性尖峰，第一极小值的位置也随之漂移。数据生成后先统一重采样或降采样再进算法。第二个坑是嵌入维数选得比实际需要大很多。Cao 方法虽然能给出一个饱和点，但超过这个点之后继续增加 m，最近邻搜索的计算量按 O(m²) 增长，而有效相点数量 M 线性减少，统计误差反而变大。第三个坑是把互信息找出来的 τ 直接套在非混沌数据上。Cao 方法的理论基础是 Takens 嵌入定理，它要求系统是确定性的，如果输入是纯随机序列，E2(d) 会一直不收敛，此时强行取某个 d 作为嵌入维数没有意义。
记住一个判断原则：当互信息曲线找不到明显的第一极小值时，先检查数据长度和数据质量，而不是急着调 bin 数或换延迟范围。Rossler 系统仿真的数据通常不会出这个问题，但如果采集的是实验信号，可能需要在数据预处理阶段先做带通滤波或去趋势，再重新计算互信息和 Cao 方法。
本文还有配套的精品资源，点击获取
