递归量化分析（RQA）的MATLAB实现：离散时间序列的递归图与七指标

简介RQA对离散时间序列进行递归图分析的MATLAB代码包面向时间序列分析与信号处理方向的学生、研究人员和工程技术人员提供完整的数据预处理、递归图构建与量化分析流程可用于心电图、脑电、金融波动和机械故障诊断等实际场景。整个zip压缩包仅1.3MB包含2个文件一个MATLAB脚本与一份PDF说明文档脚本实现了从数据读入到提取复发率、最长对角线长度、平均长度线和熵等核心指标的处理链路文档则针对每个量化指标给出计算说明与解读建议并解释距离阈值等关键参数的选择方法。目前已有160人学习下载适合作为快速入门RQA的轻量示例。通过这套代码读者可直接将其用于自己的离散序列修改参数观察指标变化快速验证非线性动态特征同时结合PDF说明理解递归图背后的理论为后续深入分析打下基础。1. 递归量化分析RQA对离散时间序列到底在算什么一张递归图和七个统计量很多从业者手里有一段离散时间序列传感器振动值、设备温度、订单量、心电图R-R间期。用均值、方差和功率谱看信号平稳得像什么都没发生可系统实际上已经换了状态。递归量化分析Recurrence Quantification Analysis的思路跟线性工具完全不同先把一维序列嵌入高维相位空间画一张递归图Recurrence Plot再用递归率、确定度、最长对角线等指标把图里的纹理压缩成数字。这份资源是一套MATLAB代码输入一列时间序列即可完成相位重构、递归图生成和RQA指标计算适合做故障诊断、状态识别或非线性信号分析的从业者也适合刚接触非线性时间序列、想拿一套能跑通的代码做对比实验的人。2. 为什么要看轨迹重访从递归图到RQA指标的构造逻辑2.1 相位空间重构和多维向量的距离一维时间序列看不出系统状态是因为单个变量只记录了系统在某一方向上的投影。Takens延迟嵌入定理告诉我们用延迟坐标构造的向量可以恢复系统的相空间拓扑结构。对序列 u(1), u(2), …, u(N)嵌入维数为 m延迟时间为 τ得到的第 i 个状态向量为x(i) [u(i), u(iτ), u(i2τ), …, u(i(m-1)τ)]注意这里的 i 取值范围是 1 到 N - (m-1)τ。这个有效点数 Nrec 是整个分析里最容易出错的地方后面避坑章会专门讲。构造好状态向量后递归图的本质就是一个阈值化的距离矩阵。对任意两个向量 x(i) 和 x(j)计算两者之间的距离如果距离小于等于临界半径 ε就认为系统在 i 时刻和 j 时刻的轨迹发生了“重访”在递归图上记 1否则记 0。距离的取法直接影响递归图纹理。最常用的是欧氏距离即 L2 范数也有人用切比雪夫距离L∞ 范数它取的是各分量差值的最大值。对同一段数据L∞ 通常会给出更大的递归区域因为它的判定条件更宽松只要有一个分量接近就算递归。我一般会把默认算法写成欧氏距离切比雪夫作为备用选项方便在不同任务里快速切换。2.2 递归矩阵与七个核心指标的物理解释递归矩阵 R 是 Nrec × Nrec 的 0/1 矩阵R(i,j) 表示 x(i) 与 x(j) 是否足够接近。对所有递归点做统计分析就得到RQA指标。最核心的七个指标分别描述递归图的不同侧面。递归率 RR 是图中递归点占所有点对的比例反映系统轨迹在相位空间中的“密度”。确定性 DET 是落在长度不小于 Lmin 的对角线段上的递归点数量占递归点总数的比例反映系统的确定性和周期性。平均对角线长度 L 是所有长度不小于 Lmin 的对角线段的平均长度最长对角线 Lmax 则是主对角线之外最长的连续对角线段。这两个指标对混沌系统中的轨道分离速度很敏感。对角线熵 ENTR 是对角线长度分布的香农熵衡量递归图纹理的复杂程度周期信号熵低混沌信号熵高。剩下两个指标面向垂直方向。层流率 LAM 是落在垂直方向连续线段上的递归点占比捕获系统“停留在某个状态附近”的时间平均垂直长度 TT 也叫湍流时间衡量系统在局部状态停留的平均时间。这七个指标组合起来相当于把一张二维纹理图压缩成一组可以跨信号比较的标量。2.3 离散序列和连续序列在递归图上的本质差异连续时间序列的距离分布是连续的可以选一个合适的 ε 让递归率落在合理范围但离散序列尤其是整数值序列距离只取少数几个值。比如取值为 {0, 1, 2, 3} 的符号序列任意两个状态向量之间的距离只有有限种可能同一个 ε 下会出现大片连续递归区域递归图看起来像一块块色斑而不是连续信号那种细碎纹理。这意味着离散序列用固定绝对阈值很容易翻车ε 设小递归点太少ε 加一递归点猛增。更稳的做法是用固定百分比邻域也就是对每个点只保留距离最近的固定比例邻域作为递归点或者用第 4 章介绍的目标递归率校准法。还有一个常见处理是先对离散序列做差分或标准化再进行相位空间重构这样能减小取值离散程度对递归图纹理的扭曲。3. MATLAB实现三段式代码把时间序列变成可比较的指标3.1 相位空间重构函数有效点数的来源代码包里第一个核心函数是相位空间重构。输入一维列向量 u返回状态向量矩阵 Y每一行对应一个相位空间点。function [Y, Nrec] phaseSpaceReconstruct(u, m, tau) % 时间延迟嵌入将一维序列映射到 m 维相位空间 % u : 一维时间序列列向量 % m : 嵌入维数正整数 % tau : 时间延迟正整数单位是采样间隔倍数 % Y : Nrec x m 的状态向量矩阵 % Nrec: 有效点数等于 N - (m-1)*tau u u(:); N length(u); Nrec N - (m-1)*tau; if Nrec 1 error(序列长度不够请减小 m 或 tau); end Y zeros(Nrec, m); for j 1:m % 第 j 列对应延迟 (j-1)*tau Y(:, j) u((j-1)*tau (1:Nrec)); end end这段代码的关键是起点和终点。每个状态向量最后一个分量是 u(i(m-1)τ)所以 i 最大只能取到 N-(m-1)τ。这里输出的 Nrec 就是后续递归矩阵的尺寸。m 和 τ 越大Nrec 越小数据浪费越多对采样点只有几百的短序列m 不建议超过 5。还要注意 u 被强制转为列向量避免行向量输入时矩阵形状出错。3.2 递归矩阵生成距离计算方式与阈值化有了状态向量矩阵下一步是构造递归矩阵。下面这个函数默认用欧氏距离同时支持切比雪夫距离。function R recurrenceMatrix(u, m, tau, eps, normType) % 生成递归矩阵 R逻辑矩阵1 表示递归点 % u : 一维时间序列 % m, tau : 嵌入维数与延迟 % eps : 临界半径标量 % normType: euclidean 或 chebychev Y phaseSpaceReconstruct(u, m, tau); Nrec size(Y, 1); R false(Nrec, Nrec); for i 1:Nrec diffMat Y - Y(i, :); switch lower(normType) case euclidean d sqrt(sum(diffMat.^2, 2)); case chebychev d max(abs(diffMat), [], 2); otherwise error(不支持的范数类型); end R(i, :) (d eps).; end end这段实现用循环逐行计算距离逻辑清晰但速度一般。对 Nrec 在 1000 以下运行没问题Nrec 到 5000 时计算时间会明显上升。如果机器上装了 MATLAB 自带的统计与机器学习工具箱可以把内层循环换成 pdist2能快一个数量级。R 矩阵用逻辑类型存储每个元素只占 1 字节这是控制内存的关键习惯不要用 double 存 0/1。3.3 RQA统计对角线和垂直线段的运行长度编码计算RQA指标的核心是对递归矩阵的每个对角线偏移和每个垂直方向做运行长度统计。下面给出一个能直接使用的指标函数function rqa computeRQA(R, Lmin) % R : 逻辑递归矩阵已去除主对角线 % Lmin : 参与统计的最小线段长度默认取 2 if nargin 2, Lmin 2; end N size(R, 1); % 递归率去掉主对角线后的点对总数 % 主对角线全部置为 0避免自递归点干扰 R R .* ~eye(N); NR sum(R(:)); totalPairs N * (N - 1); rqa.RR NR / totalPairs; % 对角线方向统计 diagFreq zeros(1, N); sumDiagLength 0; maxDiag 0; for d 1:N-1 seg R(1:N-d, 1d:N); % 这一列依次按行展开后连续为 1 的段就是一段对角线段 lineVals seg(:).; idx find(diff([0, lineVals, 0]) ~ 0); for k 1:2:length(idx)-1 len idx(k1) - idx(k); if len Lmin diagFreq(len) diagFreq(len) 1; sumDiagLength sumDiagLength len; maxDiag max(maxDiag, len); end end end rqa.L sumDiagLength / max(1, sum(diagFreq(Lmin:end))); rqa.Lmax maxDiag; rqa.DET sumDiagLength / max(1, NR); % 垂直方向统计 vertFreq zeros(1, N); sumVertLength 0; for col 1:N colVals R(:, col).; idx find(diff([0, colVals, 0]) ~ 0); for k 1:2:length(idx)-1 len idx(k1) - idx(k); if len Lmin vertFreq(len) vertFreq(len) 1; sumVertLength sumVertLength len; end end end rqa.LAM sumVertLength / max(1, NR); rqa.TT sumVertLength / max(1, sum(vertFreq(Lmin:end))); % 对角线熵基于对角线段长度分布 p diagFreq(Lmin:end) / max(1, sum(diagFreq(Lmin:end))); p p(p 0); rqa.ENTR -sum(p .* log(p)); end这段代码把每个偏移方向上的序列拉直后做运行长度统计原理上等价于沿对角线扫描。重点逻辑是 idx 变量定位连续 1 的起止位置从而实现任意长度的线段统计。DET 的分母用的是 NR也就是去除主对角线后的总递归点数如果某个实现把主对角线也算进递归点DET 会虚高跨工具对比时要注意。Lmin 默认取 2 是多数文献的惯例但如果你想复现一篇具体论文务必先确认论文中的 Lmin 取值。3.4 一次完整运行从加载数据到输出指标表和递归图把三段函数串起来主脚本只需要几行load demoSignal.mat; % 假设数据里有一个变量 sig是列向量 m 3; tau 1; eps 0.5; Lmin 2; Y phaseSpaceReconstruct(sig, m, tau); R recurrenceMatrix(sig, m, tau, eps, euclidean); R R .* ~eye(size(R, 1)); % 主动剔除主对角线 rqa computeRQA(R, Lmin); imagesc(R); title(Recurrence Plot);运行后输出 rqa.RR、rqa.DET、rqa.L、rqa.Lmax、rqa.ENTR、rqa.LAM、rqa.TT。我习惯把 R 矩阵直接显示成黑白图先肉眼看纹理再做数值统计因为很多异常信号在递归图上会有明显特征突变表现为横向或纵向的断续线漂移表现为递归图边缘变暗。建议每一步都打印 Nrec确认有效点数和预期一致。4. 三个关键参数怎么设嵌入维数、时间延迟与临界半径的实操4.1 嵌入维数m用Cao方法代替人工试错嵌入维数 m 决定状态向量的长度。m 太小相空间没有充分展开距离计算掩盖系统结构m 太大每个状态向量需要的延迟点数变多Nrec 变小统计噪声变大。实际工作里很少直接试 m我一般用 Cao 方法它不需要先验地估计噪声水平适合短序列。Cao 方法的核心是计算 E1(d) E(d1) / E(d)其中 E(d) 是基于最近邻距离定义的均值。当 E1(d) 随着 d 增大开始逼近 1 并保持平稳时对应的 d 就是合适的嵌入维数。下面的函数计算 E1 随维数的变化function E1 caoE1(u, maxDim, tau) % 简化版 Cao 方法只输出 E1 指标 u u(:); for d 1:maxDim Y phaseSpaceReconstruct(u, d, tau); Nrec size(Y, 1); % 对每个点找最近邻距离 nnDist zeros(Nrec, 1); for i 1:Nrec dst sqrt(sum((Y - Y(i,:)).^2, 2)); dst(i) inf; % 排除自身 nnDist(i) min(dst); end a zeros(Nrec, 1); % 计算加入下一维后的距离变化 Y1 phaseSpaceReconstruct(u, d1, tau); if size(Y1, 1) Nrec Nrec size(Y1, 1); Y Y(1:Nrec, :); nnDist nnDist(1:Nrec); end for i 1:Nrec dstPlus sqrt(sum((Y1 - Y1(i,:)).^2, 2)); dstPlus(i) inf; nnPlus min(dstPlus); a(i) dstPlus / max(nnDist(i), eps); end E(d) mean(a); end E1 E(2:end) ./ E(1:end-1); plot(2:maxDim, E1, -o);这段代码没有做文献里的分段平均改进但对大多数中等长度序列足够用。实际操作时我会把 E1 从 2 画到 8观察曲线第一次进入“平台区”的转折点。例如 E1 在 d4 之后稳定在 1 附近且没有继续明显上升就选 m4。Cao 方法输出的不是一个唯一确定值而是一个合理区间所以不必纠结分数差异。4.2 时间延迟τ互信息极小值与自相关法τ 决定每个状态向量两个相邻分量的间隔。τ 太大相邻分量几乎不相关重构轨迹被拉散τ 太小分量之间冗余递归图纹理被压缩。最常用的两套选择法自相关法和互信息法。自相关法简单直接计算序列自相关函数找到它首次下降到初始值的 1-1/e 约 36.8% 时的滞后作为 τ。互信息法更全面它计算原始序列和延迟序列之间的互信息 I(τ)取第一个极小值对应的 τ。互信息能捕捉非线性相关性所以对混沌信号更可靠但计算量稍大。一个简化版的互信息计算如下function I mutualInfo(u, tau) % 简化版互信息用直方图估计联合分布 u u(:); N length(u); x u(1:N-tau); y u(1tau:N); nbins max(5, round(sqrt(numel(x)))); [c, edges] histcounts(x, nbins); px c / sum(c); [~, edgesY] histcounts(y, nbins); py histcounts(y, edgesY) / numel(y); pJoint histcounts2(x, y, edges, edgesY) / numel(x); logTerm pJoint ./ (px(:) * py(:).); logTerm(pJoint 0) 1; % 避免 log(0) I nansum(pJoint .* log(abs(logTerm))); end这个实现没有做互信息的均值校正但找极小值的趋势是可靠的。我的习惯是把 τ 从 1 到 10 都算一遍画出 I(τ) 曲线取第一个明确极小值。如果曲线太平、没有清晰极小值就用 τ1 作为默认值因为对采样率较高的离散信号过大的 τ 会丢失时间结构。4.3 临界半径ε固定阈值、百分比阈值与目标递归率ε 是 RQA 里最敏感的参数它直接决定递归点是稀疏还是密集也影响 DET 和 Lmax。经验上ε 太小递归点太少对角线统计噪声大ε 太大递归点太多几乎全部点都变成递归点DET 趋近 1 失去区分度。说直白点ε 的选取在有些论文里写得像玄学实质是尺度选择问题。三种常用方案各有边界。固定绝对阈值适合所有被比较的序列具有相同量纲和幅值的场景比如同一传感器的多段振动数据百分比阈值是把每个点到其所有近邻的距离排序取某个分位数作为该点的自适应阈值适合幅值漂移信号目标递归率校准是先设定期望的 RR比如 0.05再用二分法搜索 ε使实际 RR 逼近目标值。第三种方案对跨工况对比最友好我会在最后一章给出具体代码。4.4 常用参数组合表与默认建议下面这个表是跨多个使用场景的经验区间不是硬性标准用于给第一次跑这套代码的人一个参考起点数据形态嵌入维数 m延迟 τε 策略连续传感器振动信号3~51~3可用互信息目标 RR 0.05~0.10离散整数计数序列3~51百分比阈值 5%~10%金融收益率序列2~31固定百分位邻域生理信号心电、脑电4~6视采样率而定目标 RR 0.02~0.08拿到新数据时我会先用 m4、tau1、目标RR0.05 作为默认粗跑一遍再根据递归图纹理和具体任务调整。不要盲目追求某个“最优”参数RQA 指标在 m 和 τ 的合理邻域内变化通常不大真正会导致指标剧烈变化的是 ε所以参数优先级是 ε τ m。5. 避坑指南在离散序列上跑RQA的六个翻车现场5.1 主对角线剔除与否影响RR和DET现象同一段数据用两套代码算 RR 和 DET结果差别很大一套 RR0.04另一套 RR0.06。原因主对角线上 x(i) 和它自身的距离恒为 0只要 ε0 就是递归点。如果实现把主对角线计入递归点总数RR 和 DET 的分母分子同时变化指标虚高不同实现是否剔除主对角线不一致导致跨工具对比直接失真。解决在统计前强制把主对角线置 0也就是代码中的R R .* ~eye(N)并且把分母固定为 N(N-1)。这个决定需要在实验记录里写明避免自己和读者误判。5.2 序列长度不同指标能不能比现象A 序列 1000 点、B 序列 500 点B 的 DET 明显更高看起来 B 更“确定”。原因短序列的有效点数少边界效应影响大递归图边缘部分不完整同时 Lmax 天然受序列长度限制短序列的 Lmax 相对噪声更大。解决比较前把各序列截取成相同长度或采用滑动窗口对全序列提取指标曲线。如果必须在不同长度下比较优先使用 RR、DET 和 LAM 这些比值型指标Lmax 只能作为辅助且最好用 Lmax/Nrec 归一化。5.3 大N递归矩阵的内存瓶颈现象Nrec 超过 5000 时计算变慢Nrec 超过 8000 时 MATLAB 内存占用巨大甚至直接无响应。原因递归矩阵是 Nrec × Nrec 的逻辑矩阵Nrec8000 时也有 64M 个元素如果中途用 double 保存距离矩阵内存占用立刻翻 8 倍。解决只保存 R 的逻辑矩阵不保留完整距离矩阵距离计算用分块或逐行处理。如果只看指标不需要绘图可以逐条偏移计算对角线统计跳过整矩阵存储。代码包里专门提供了一个内存友好版本用分块循环替代全矩阵构造。5.4 离散整数序列的“块状递归图”现象取值为 {0, 1, 2} 的序列画递归图大范围区域整块变黑几乎看不出纹理。原因离散取值少状态向量重复率高用绝对距离阈值时很多不同状态的距离也小于 ε递归点大面积连成片。解决先做序列标准化或一阶差分得到连续分布的数据再用百分比阈值。对符号序列要明确说明结果属于“符号化递归”与连续信号的 RQA 指标不在同一尺度上不能直接混比。5.5 DET高不等于混沌周期与混沌的区分现象纯正弦信号的 DET 接近 1Logistic 混沌映射的 DET 也能到 0.9 以上单看 DET 完全分不开。原因DET 只反映轨迹返回的确定性结构周期性数据和混沌数据都能产生明显的对角结构。它们真正的差异在 Lmax 和 ENT R周期信号最长对角线接近 Nrec熵低混沌信号最长对角线较短对角线长度分布更随机熵较高。解决比较时同时看 Lmax/Nrec 和 ENTR。周期信号两者呈“长而低熵”组合混沌信号呈“短而高熵”组合。如果只报 DET结论说服力不够。5.6 工具箱依赖导致代码跑不起来现象代码下载后直接运行报错提示缺少某个函数或工具箱。原因部分公开实现依赖统计工具箱里的 pdist2、图像处理工具箱里的连通域函数而基础版 MATLAB 不带这些模块。解决资源内的实现全部基于基础 MATLAB不依赖任何附加工具箱。如果要用 pdist2 加速运行前先判断license(test, Statistics_Toolbox)不可用时回退到循环版本。这样既保证可复现又能照顾有工具箱的用户。6. 进阶技巧滑动窗口RQA与目标递归率校准6.1 滑动窗口RQA的快速实现RQA 指标的静态值只适合平稳段分析对非平稳信号和状态切换检测需要把序列切成窗口滑动计算。窗口长度建议取 200 到 1000 点步长取窗口的 1/4 到 1/2。核心思路是对每个窗口重新做相位空间重构和 ε 校准得到一组随时间变化的指标序列。function rqaSeries slidingRQA(u, W, S, m, tau, targetRR) tAxis W:S:length(u); rqaSeries zeros(numel(tAxis), 7); k 0; for t tAxis k k 1; seg u(t-W1:t); epsOpt calibrateEpsilon(seg, m, tau, targetRR); R recurrenceMatrix(seg, m, tau, epsOpt, euclidean); R R .* ~eye(size(R, 1)); rqa computeRQA(R, 2); rqaSeries(k, :) [rqa.RR, rqa.DET, rqa.L, rqa.Lmax, rqa.ENTR, rqa.LAM, rqa.TT]; end end6.2 校准ε让不同信号在同样递归密度下比较对每个窗口重新校准 ε是为了保证滑动过程中 RR 始终接近目标值这样才能把 DET、Lmax 的变化归因于动力学结构变化而不是递归密度变化。二分法搜索在 ε 与 RR 单调递增的前提下可靠迭代 20 次足够收敛。function epsOpt calibrateEpsilon(u, m, tau, targetRR) Y phaseSpaceReconstruct(u, m, tau); Nrec size(Y, 1); % 计算距离上界用于设置二分区间 maxDist 0; for i 1:Nrec dst sqrt(sum((Y - Y(i,:)).^2, 2)); maxDist max(maxDist, max(dst)); end lo 0; hi maxDist; for iter 1:20 mid (lo hi) / 2; R recurrenceMatrix(u, m, tau, mid, euclidean); RR sum(R(:) ~eye(Nrec)) / (Nrec * (Nrec - 1)); if RR targetRR lo mid; else hi mid; end end epsOpt (lo hi) / 2; end这段代码会随递归矩阵计算重复运行代价不小但胜在稳定。我之前用固定 ε 对比两段不同工况的振动数据结果一个 RR 是 0.002另一个是 0.15DET 完全不可比后来写成目标递归率校准才发现两段信号在动力学复杂度上其实很接近。从那以后凡是需要跨工况或跨样本比较 RQA 指标我都强制走一遍“定窗口长度—校准 RR—检查递归图边缘是否有截断”的流程希望帮到你。本文还有配套的精品资源点击获取