基于Rosenstein小数据量法的最大Lyapunov指数计算源码

资 源 简 介

本项目利用MATLAB编程语言实现Rosenstein算法（亦称小数据量法），专门用于从短时间序列或含有噪声的实验数据中计算最大Lyapunov指数（LLE），以判定系统的混沌特性。核心功能流程包括：首先对输入的单变量时间序列进行预处理（去噪、去趋势、归一化）；其次，依据Takens嵌入定理进行相空间重构，程序内置了辅助算法（如互信息法计算延迟时间tau，虚假邻点法计算嵌入维数m，或集成C-C法）来自动或交互式确定最优重构参数；接着，利用快速傅里叶变换（FFT）计算时间序列的平均周期，用于在寻找最近邻点时排除时间相关性带来的同轨误差；算法核心部分将在重构的高维相空间中能够为每个参考向量寻找其欧氏距离最小的最近邻点，并计算这两点随着时间步长增加后的距离发散情况。系统最终输出平均对数发散度$<ln(dj)>$随时间演化的曲线图，并允许用户通过图形界面或自动算法选择曲线中的线性增长区域（线性标度区），利用最小二乘法拟合该区域的斜率，从而得到最大Lyapunov指数的估计值。该项目适用于样本数据量较小（如少于5000点）的非线性动力学分析场景。

详 情 说 明

基于Rosenstein小数据量法的最大Lyapunov指数计算系统

项目简介

本项目是一个基于MATLAB开发的非线性动力学分析工具，专门用于计算时间序列的最大Lyapunov指数（Largest Lyapunov Exponent, LLE）。项目采用经典的 Rosenstein算法（小数据量法），特别适用于数据量较少、包含噪声的实验数据或混沌系统仿真数据。

系统集成了从数据生成、预处理、相空间重构参数自适应估计、核心算法计算到结果自动拟合与可视化的全流程。代码通过分析相空间中最近邻点对随时间的距离发散率，量化系统的混沌特性。

功能特性

• 混沌数据生成：内置Lorenz混沌系统数值模拟功能，使用四阶Runge-Kutta法（RK4）生成标准测试数据。
• 数据预处理：自动剔除瞬态数据，并执行Z-score标准化（去均值、单位方差），提高计算稳定性。
• 参数自适应估计：

* 平均周期计算：利用快速傅里叶变换（FFT）计算信号平均周期，用于确定Theiler窗大小，消除去时间相关性。 * 延迟时间 ($tau$)：基于互信息法（Mutual Information），自动搜索互信息函数的第一个极小值点。 * 嵌入维数 ($m$)：基于虚假邻点法（False Nearest Neighbors, FNN），寻找虚假邻点比例低于阈值（5%）的最小维数。

• Rosenstein核心算法：

* 高效的相空间重构与最近邻搜索。 * 引入Theiler窗机制，避免时间序列自身相关性造成的同轨误差。 * 计算并累积平均对数发散度 $langle ln(d_j) rangle$。

• 自动线性拟合：采用滑动窗口搜索策略，自动识别平均对数发散曲线中的最佳线性标度区（基于最大 $R^2$ 原则），利用最小二乘法估算LLE。
• 多维可视化：生成综合分析图表，包括原始时序图、参数估计指标图、三维相空间轨迹图以及对数发散曲线拟合图。

系统要求

• MATLAB R2016b 或更高版本。
• 不需要额外的工具箱（代码主要依赖MATLAB基础函数，findpeaks 等函数通常包含在信号处理工具箱中，但基础数学运算均为原生实现）。

使用方法

1. 确保代码文件在MATLAB路径中。
2. 直接运行 main 函数。
3. 系统将在命令窗口输出计算过程中的状态信息（平均周期、最优 $tau$、最优 $m$、最终LLE估计值等）。
4. 运行结束后将弹出一个包含四个子图的图形窗口，展示分析结果。

---

核心代码逻辑与实现细节

本项目的主要逻辑在 main.m 中实现，具体流程如下：

1. 数据准备与预处理

程序首先通过 ode_rk4 函数求解Lorenz方程组。设定参数 $sigma=10, rho=28, beta=8/3$，积分步长 dt=0.01s。为了消除初始条件的影响，程序自动丢弃生成的轨迹的前2000个数据点。随后，通过减去均值并除以标准差的方式对数据进行归一化处理。

2. 平均周期的计算 (Theiler窗)

为了在寻找最近邻点时排除时间上相邻点的干扰，程序调用 get_mean_period_fft 子函数。该函数对去均值后的数据进行FFT变换，找到频谱幅值最大的频率（忽略直流分量），其倒数即为平均周期 mean_period。此值随后被用作Rosenstein算法中的Theiler窗参数 w。

3. 重构参数的自动估计

• 延迟时间 (tau)：程序计算数据在不同滞后下的互信息值，并通过 findpeaks 函数寻找互信息曲线取负后的峰值，从而定位第一个极小值点作为最优延迟时间。如果未找到极小值，则使用默认值10。
• 嵌入维数 (m)：程序计算不同维数下的虚假邻点比例（FNN Ratio）。通过逻辑判断找到第一个FNN比例小于0.05（5%）的维数索引，并据此确定最优嵌入维数。对于Lorenz系统，理论上通常得到 $m approx 3$。

4. 相空间重构与距离演化计算

这是Rosenstein算法的核心部分：

• 构建矩阵 $Y$（大小为 $M times m$）作为重构的相空间。
• 遍历相空间中的每一个参考点 $j$。
• 最近邻搜索：在寻找最近邻点时，严格执行Theiler窗约束，即排除索引在 $[j-w, j+w]$ 范围内的点，仅在有效范围内搜索欧氏距离最小的邻点 nearest_idx。
• 发散度追踪：计算参考点与最近邻点在后续 max_steps 个时间步长内的距离演化。为了计算效率，初始搜索使用距离平方，后续演化计算实际距离。
• 对数平均：将所有参考点对在相同演化步长 $k$ 下的距离取对数并累加，最终求得平均对数发散度序列 avg_log_dist。

5. 线性区域自动识别与LLE计算

程序不依赖人工选取线性区，而是实现了一个自动化算法：

• 使用双重循环遍历对数发散曲线的不同起始和结束位置（滑动窗口）。
• 对每个窗口内的数据进行一次多项式拟合（线性回归）。
• 计算拟合的决定系数 $R^2$。
• 选择 $R^2$ 最高且斜率为正的窗口作为“最佳线性标度区”。
• 该最佳区域的斜率即为输出的 最大Lyapunov指数 (LLE)。

6. 结果可视化

程序最后绘制一个 2x2 的图表：

• 左上：预处理后的时间序列波形。
• 右上：互信息随延迟变化的曲线（标注最优 $tau$）以及虚假邻点比例随维数变化的曲线。
• 左下：重构相空间的轨迹图（如果是3维及以上，绘制3D图；否则绘制2D图）。
• 右下：平均对数发散度随演化时间的曲线，并在图上用红色虚线标出自动拟合的线性区域，标题中显示计算得到的LLE值。

---

关键子函数说明

以下函数在代码内部定义并支撑主流程：

• ode_rk4：实现了标准的四阶Runge-Kutta数值积分算法，用于求解常微分方程组。
• lorenz_eq：定义了Lorenz系统的动力学方程，即 $dot{x}, dot{y}, dot{z}$ 的导数关系。
• get_mean_period_fft：执行FFT变换，分析频谱主峰，将频率转换为以“数据点数”为单位的平均周期。
• compute_mutual_information：计算时间序列与其延迟序列之间的互信息，利用直方图法估计概率分布。（注：代码片段显示该函数调用存在，用于tau的估计）。