简介面向时间序列分析与信号处理研究者的MATLAB递归量化分析RQA工具包用于对离散时间序列构建递归图并提取复发率、最长对角线、平均线段长度和熵等量化指标适合需要考察系统非线性动态、周期性或状态切换的科研与工程场景。压缩包共2个文件包含一个可直接运行的.m脚本与一份PDF说明文档整体仅1.3MB轻量易用。已有160人学习下载。通过示例代码可快速了解数据预处理、相似性阈值设定、递归图构建到特征量化的完整流程并参考说明文档解读各指标含义从而迁移到心电图、脑电、金融波动或机械振动等实际数据上。对于具备基础MATLAB操作能力、希望掌握非线性时间序列分析方法的从业者这份代码包能提供清晰可复现的入门参考与二次开发起点。1. RQA 递归量化分析不是玄学离散时间序列能挖出多少隐藏结构某次数据分析任务里一条来自振动台架的传感器信号频谱上几乎看不出故障前后的差异功率谱密度曲线看起来都很平坦我尝试用递归量化分析RQA把离散时间序列重构成高维相空间计算递归图中对角线结构的比例故障前后 DET 指标差了近一倍问题才一下子暴露。这类代码包在 MATLAB 生态里流传很广标题里的“RQA对离散时间序列进行递归图分析.zip”就是其中一种常见实现。RQA 的核心不是画一张花哨的递归图而是把“系统轨迹有没有回到相似状态”这件事量化为一串指标适合做故障诊断、状态识别、混沌特征判断。本文按“先懂原理、再抄代码、后避坑”的顺序带你把它真正落地而不是停留在跑通 demo 的阶段。2. 递归图与 RQA 指标先搞懂 D 维空间里的“点与线”2.1 从离散序列到相空间重构嵌入维 m 和延迟 τ 怎么定离散时间序列本质上是一串标量观测值比如每隔 1 毫秒采一次电压、每个交易日记一个价格。单独看这些标量只能看到幅度变化看不出系统内部有多少自由度。RQA 的第一步就是把这串标量放到相空间里去做法是构造延迟向量X(i) [x(i), x(iτ), x(i2τ), …, x(i(m-1)τ)]其中 m 是嵌入维τ 是延迟步数i 取值范围是 1 到 N-(m-1)τ。重构后每个时刻得到一个 m 维点这些点在相空间里连成一条轨迹递归图记录的就是这条轨迹什么时候回到某个点的邻域。这个思想来自时间延迟嵌入定理实际操作中没有人会去严格验证维数而是靠指标选 m 和 τ。我一般先算平均互信息AMI随 τ 的变化取第一个局部极小值作为延迟再用假近邻FNN或 Cao 方法定 m。MATLAB 里没有内置的 AMI 函数但实现不复杂把序列分成若干等宽区间统计联合概率和边际概率计算互信息。还有更省事的替代看自相关函数取自相关第一次衰减到初始值 1/e 附近的延迟。这个做法对过采样的信号尤其有效因为过密采样会让相邻样本高度相关τ1 时递归图全是“厚对角线”好像确定性很强其实是采样率制造的假象。血泪经验是不要为了“包含更多历史信息”而把 m 设得很大。有一段 N200 的短序列我一开始用 m10重构后向量只剩 110 个距离矩阵稀疏到没法看DET 和 LAM 几乎随机跳动。后来把 m 压回 3指标才稳定。离散时间序列长度有限时m 越大有效向量越少不确定性越高我自己在工程里常用范围是 m2 到 7短序列建议从 m2、τ1 起步先看递归图形态再微调。Cao 方法的好处是它不依赖阈值比 FNN 稳定适合批处理。2.2 RQA 核心指标一览RR、DET、L、LAM、ENTR 到底在说什么有了递归矩阵 RR(i,j)1 表示第 i 个重构向量和第 j 个重构向量之间的距离小于阈值 ε。矩阵里全是点和空白人眼能看出纹理但机器要的是数字。RQA 指标就是把纹理压缩成一组统计量下面这几个是我最常用的。指标全称/计算对象高值说明RR递归点密度R 中 1 的占比状态重现的总频率过高多半是阈值太大DET大于等于 lmin 的对角线段上的点占全部递归点的比例轨迹遵循相似路径的比例代表确定性结构L所有对角线段长度的均值系统能保持相似演化路径的平均长度Lmax最长对角线段长度存在周期或很稳定的结构但可能是单根离群段LAM大于等于 vmin 的垂直线段上的点占全部递归点的比例轨迹短暂“卡”在某个区域的比例反映间歇与层流TT垂直线段平均长度系统停留在同一状态的平均时间ENTR对角线段长度分布的 Shannon 熵递归结构越复杂熵越高计算 DET 最常见的公式是DET Σ_{l≥lmin} l·P(l) / Σ_{l≥1} l·P(l)其中 P(l) 是长度为 l 的对角线段出现次数。LAM 的定义类似只不过把对角线换成垂直线。从公式能看出lmin 和 vmin 直接决定了指标口径通常在 2 到 5 之间取整数。对量化很糙的离散序列比如只取 0 和 1 的符号序列长度 1 的线段会特别多DET 会被压得很低这时候不要立刻怀疑代码先确认最小线长是否合理。需要特别提醒RQA 指标对离散序列长度非常敏感。N 在 1000 以上时 RR、DET 的估计才比较稳N 在 200 到 500 时单次计算结果波动大最好用滑动窗口或多次重复采样取平均。不要拿着 100 个点的序列去断言“系统是混沌的”那是在赌运气。2.3 阈值 ε 的选择固定阈值、百分比阈值与最大距离阈值RQA 里最容易被说成“玄学”的就是阈值 ε。“为什么你取 0.2不取 0.15”这个问题我也被问过很多次。ε 本质上是在定义“多近算相似”它必须与距离的量纲和数据动态范围一致。三种常见选法各有利弊。固定绝对阈值实现最简单但一遇到量纲不同的数据就翻车。序列 A 振幅是 0.01 量级序列 B 是 100 量级用同一个 ε 根本没有可比性。所以固定阈值之前我几乎总是先把序列标准化为 z-scorex (x - mean(x)) / std(x)。标准化之后ε 取 0.1 到 0.3 通常能给出合理密度但也要看数据形状。百分比阈值是按距离分布取分位数。比如想要递归点密度大约 15%就先算所有重构向量两两之间的距离取第 85 百分位距离作为 ε。这种做法在不同数据集之间可比性最好因为每个数据集的 RR 都被拉到同一水平。缺点是它掩盖了“原本状态重现就不频繁”的差异因此用于比较系统复杂度时没问题用于异常检测时要谨慎。最大距离比例阈值是把 ε 定义为 α 乘以距离矩阵最大值α 常见 0.05 到 0.3。这个方案不受量纲影响但最大值容易被个别离群距离带偏。我的习惯是先把距离矩阵拍平画直方图看低尾位置再综合选择。小技巧是用下面这段快速找出分位数阈值d D(:); d sort(d); idx round(0.15 * numel(d)); % 期望约15%的递归点 eps d(max(idx, 1));这里的 d 是距离矩阵展开后的向量D 来自第 3 章的递归矩阵函数。用分位数阈值有一个附带好处递归点密度不再依赖数据幅值报告结果时也更容易复现。无论选哪种请一定把 ε 和选法写进文档否则三个月后你连自己都看不懂自己算的是什么。3. MATLAB 代码落地把 RQA 跑通的最小实现与参数设置标题里的代码包解压后我认为没必要逐行去读里面的每个画图脚本先把两个核心函数理清楚一个算递归矩阵一个算 RQA 指标。我习惯把这两件事拆开因为调参时只需要换 ε、m、τ不需要碰指标计算。下面这套骨架是我在类似任务里反复用的你可以直接对照着重构自己的版本。3.1 用 MATLAB 从零写递归矩阵最小函数与内存说明function R recurrenceMatrix(x, m, tau, eps) % x: 输入离散时间序列一维列向量 % m: 嵌入维数正整数 % tau: 延迟步数正整数 % eps: 距离阈值建议在标准化后的尺度下选择 x x(:); x (x - mean(x)) / std(x); % 先做 z-score 标准化 N length(x); M N - (m - 1) * tau; % 重构向量个数 if M 2 error(数据太短重构后少于2个向量); end % 相空间重构X 的每一行是一个 m 维向量 X zeros(M, m); for i 1:M X(i, :) x(i:tau:i (m - 1) * tau).; end % 两两距离的平方利用 |a-b|^2 |a|^2 |b|^2 - 2ab X2 sum(X.^2, 2); D2 X2 - 2 * (X * X.) X2.; D2(D2 0) 0; % 消除浮点运算引入的微小负值 D sqrt(D2); R D eps; % 距离小于阈值记为递归点 R(1:M1:end) false; % 去掉主对角线自匹配 end这段代码的思路是先用矩阵运算一次性算出所有状态向量的两两欧氏距离再按阈值比较。标准化这一步不是可选的除非你明确知道数据量纲一致且要做固定阈值比较。X 的行数是 M列数是 mX * X. 是 M×M 的矩阵乘法当 M 很大时会占内存M2000 时距离矩阵约 32 MBM10000 时直接到 800 MB。所以长序列请务必分段滑动或者把距离计算改成循环逐块计算不要在一个函数里硬撑。关于阈值参数 eps注意它现在是标准化后 z-score 尺度下的绝对距离。如果前面做了标准化eps0.2 表示两个状态向量的欧氏距离小于 0.2 倍标准差。如果你更想用百分比阈值可以在调用这个函数之前算好 eps再把 eps 传进来。参数 m 和 tau 的选择见 2.1不要再让主程序里出现魔法数字。3.2 计算 RQA 指标从递归矩阵到 RR、DET、LAM、熵function out rqaMetrics(R, minLen) % R: 递归矩阵逻辑矩阵或0/1矩阵 % minLen: 最小对角/垂直线长度默认2 if nargin 2, minLen 2; end N size(R, 1); R double(R); out.RR sum(R(:)) / numel(R); % 递归点密度 % ---- 对角线段长度统计 ---- diagCounts zeros(N, 1); for d -(N-1):(N-1) seg diag(R, d) 1; L runLengths(seg); for k 1:length(L) diagCounts(L(k)) diagCounts(L(k)) 1; end end totalDiagPoints sum(diagCounts .* (1:N).); out.DET sum(diagCounts(minLen:N) .* (minLen:N).) / totalDiagPoints; out.L totalDiagPoints / max(sum(diagCounts), eps); idx find(diagCounts 0); if isempty(idx) out.Lmax 0; else out.Lmax idx(end); end p diagCounts / max(sum(diagCounts), eps); p p(p 0); out.ENTR -sum(p .* log2(p)); % 对角长度分布的信息熵 % ---- 垂直线段长度统计 ---- vertCounts zeros(N, 1); for j 1:N L runLengths(R(:, j) 1); for k 1:length(L) vertCounts(L(k)) vertCounts(L(k)) 1; end end totalVertPoints sum(vertCounts .* (1:N).); out.LAM sum(vertCounts(minLen:N) .* (minLen:N).) / totalVertPoints; out.TT totalVertPoints / max(sum(vertCounts), eps); end function L runLengths(seg) % 统计一个逻辑向量中连续1段的长度 seg seg(:); b [false; seg; false]; startIdx find(diff(b) 1); endIdx find(diff(b) -1); L endIdx - startIdx; end这个函数把递归矩阵拆成对角线和垂直线两大类结构。ripple 在于 runLengths 统计的是每个最长连续段的长度例如某条对角线上有连续 5 个点就计一条长度为 5 的线段不会重复拆成长度 1 到 5 的多段。因此 DET 的分母 totalDiagPoints 等于所有落在线段上的递归点数分子只保留长度不小于 minLen 的段上的点数两者比值不会大于 1。使用时注意 minLen 的语义它同时作用于对角线和垂直线。标准 RQA 工具里lmin 和 vmin 可以分别设置我这里为了简洁共用一个参数如果数据里孤立点特别多建议把 minLen 设成 3 再看 DET否则长度 1 的孤立点会把 DET 拉低。输出 out 是一个 struct包含 RR、DET、L、Lmax、ENTR、LAM、TT 七个字段够大多数场景使用。3.3 参数组合表不同场景下的推荐起点没有任何一组参数能通吃所有离散时间序列但下面这张表可以作为调参起点。所有场景都默认先做 z-score 标准化。数据特征mτεminLen备注长连续信号N10002~4AMI 第一极小值0.1~0.32先看递归图纹理再微调强噪声背景3~5偏小避免过采样百分比 20%3提高 DET 抗噪能力短序列N≈2002~31~3百分比 15%~30%2不要追求高精度指标明显周期信号2~4不大于主周期一半相对宽松2DET 会非常高正常现象0/1 符号序列1~210等值递归2欧氏距离不适合可用 Hamming 距离自测代码我建议用 Logistic 映射生成一段混沌序列验证整个流程是否能得到合理结果。只要 DET 明显高于同长度白噪声的 DET就说明代码链路是通的。rng(1); r 3.8; x0 0.4; N 1000; x zeros(N, 1); x(1) x0; for i 2:N x(i) r * x(i-1) * (1 - x(i-1)); % 一维混沌映射 end m 2; tau 1; eps 0.25; minLen 2; R recurrenceMatrix(x, m, tau, eps); out rqaMetrics(R, minLen); disp(out);如果 out.DET 非常低先检查阈值是否太严递归矩阵可能几乎全空如果 DET 接近 1又被对角线丛覆盖则多半是 m 过大或阈值过大。混沌映射的 DET 通常显著高于白噪声但不会接近 1这样一组自测能帮你快速定位是重构问题还是指标计算问题。4. 离散序列 RQA 避坑指南重复值、嵌入维与阈值选错的补救办法4.1 重复值导致 RR 虚高到离谱现象传感器数据经过四舍五入后出现大量重复值比如某一档位电压在序列里连续出现几十次计算出的递归图是一块块矩形色块RR 冲到 50% 以上看起来系统“天天回到同一个状态”但这其实是量化精度造成的假象。原因绝对距离阈值下完全相等的两个点距离是 0无论 ε 多小它们都算递归点。离散时间序列如果来源是数值仿真或 ADC 采样量化台阶会制造系统性重复。解决最稳妥的办法是在标准化之前检查 x 的唯一值数量如果重复值比例超过 5%要么提高数据精度要么给信号加一个幅度很小的抖动比如随机噪声振幅取数据标准差的 1e-6。抖动不会改变动力学结构但能打破完全相等的假递归。这个方法有“后悔药”的属性我一般在最后报告里注明“已加微幅扰动”避免别人质疑。4.2 嵌入维 m 过大导致递归图变成“对角线丛”现象把 m 从 2 一直调到 6递归图越来越密DET 逼近 0.95对角线又粗又长内心觉得“终于找到强确定性了”但换一段白噪声也能得到类似结果。原因m 增大后重构向量个数 MN-(m-1)τ 快速下降。向量少了距离矩阵的基数变小只要阈值稍微给大一点任意两个向量都容易“相似”于是出现伪对角线结构。解决把 M 纳入报告。一个粗糙经验是 M 至少要大于 50且 m 不要超过 N/(3τ)。更正规的是先用 Cao 方法确定 mCao 方法会输出 E1(m) 和 E2(m) 两条曲线E2 从接近 1 变成明显大于 1 的位置就是建议的 m。很多 MATLAB 实现是把 Cao 和递归图分开算我也是这么做。4.3 阈值 ε 选择不当算出来的指标全是黑匣子现象同一段数据ε 从 0.05 调到 0.5RR 从 3% 涨到 80%DET 从 0.2 涨到 0.9。把两次结果放在同一张图里结论完全相反于是有人开始把 RQA 称为“玄学”。原因ε 本质上是判断“相似”的尺度不同尺度下递归结构不同。这不是 RQA 的缺陷而是使用者没有把阈值固定下来。解决不要只算一次我习惯做 ε 扫描设置一组候选阈值计算对应的 RR 和 DET画成曲线。曲线出现明显拐点的区域通常是一个相对稳定的 ε 段如果曲线陡峭说明数据本身对尺度极敏感这时候报告结果就要特别说明。另一个可行方案是使用百分比阈值把 RR 固定住不同数据集之间就能在统一密度下比较结构差异。无论哪种方案都必须写清楚阈值策略。4.4 白噪声和周期信号都可能有较高的 DET别急着下结论现象对一段白噪声做 RQADET 也有可能到 0.1 甚至 0.2对一段理想正弦波DET 直接接近 1。于是有人看到 DET 高就声称“有确定性”翻车概率很大。原因白噪声的随机游走也能偶相遇成长度为 2 或 3 的对角线段在递归点总数不多时占比不小正弦波本身就是周期信号轨迹无限重复DET 高是必然但不能说明它“混沌”或者“具有非线性特征”。解决加一层替代数据检验。把原始序列随机打乱生成 100 组替代序列分别计算 DET看原始 DET 是否超过 95% 分位。如果连打乱后都差不多说明统计量不敏感如果明显高于才有底气说“该序列的递归结构不是偶然”。替代数据方法不完美但它能拦住最明显的过度解读。4.5 边界效应与主对角线处理方式影响指标现象递归图四个角边缘纹理明显异常主对角线附近总是有一条特别连续的亮线把主对角线置 0 之后原本长度 3 的线段被切断成两段长度 1DET 骤降。原因重构向量来自序列前部和后部时包含的样本点数量不完整主对角线表示同一个状态与自身的距离为 0如果不去除每个点都自递归也会强化对角线结构。解决报告里写明“主对角线已剔除”或“未剔除”。剔除主对角线后最小线长 minLen 通常需要设 2 或 3否则原本跨越主对角线的连续段会被误伤。处理边界效应时我一般会丢弃首尾各 (m-1)τ 个样本对应的重构向量或者至少在结果中说明计算的 M 值。很多现成工具默认剔除主对角线但如果你用的是另一套实现务必先确认它的行为。5. 让递归图说人话可视化、置换检验与参数报告规范5.1 用 MATLAB 画递归图imagesc 与 spy 的选择递归图本身是 M×M 的 0/1 矩阵最直接的可视化是用 imagesc 把矩阵画成图像。我这里会在画图前把坐标轴方向和比例尺调好否则图像上下颠倒影响阅读。figure; imagesc(1:M, 1:M, R); colormap(gray); axis xy; axis square; xlabel(状态索引 i); ylabel(状态索引 j);M 是重构向量个数。axis xy 让 y 轴从下往上递增和矩阵显示习惯一致axis square 保证图像是正方形避免递归图纹理被拉伸变形。当 M 超过 1500 时imagesc 会明显变慢这时可以用 spy 画散点对显示速度提升很大。figure; spy(R, k, 1); axis square; title(Recurrence Plot (spy view));spy 的缺点是看不出递归密度灰度只显示点位置但它适合快速浏览大矩阵。对离散时间序列来说递归图的白色区域和黑色纹理同样重要所以不要把图缩得太小否则许多细节显示不出来。保存图片时建议输出为矢量图格式因为递归图点阵缩小时会出现严重摩尔纹。5.2 用置换检验给 DET 一个可信区间RQA 指标单独看没有绝对阈值直接说“DET0.62 表示确定性很强”是不负责的。我习惯对同一段序列做一次替代数据检验用随机打乱的时间序列作为零分布。这个过程看起来简单但能挡掉很多伪结论。nPerm 100; baseline zeros(nPerm, 1); for b 1:nPerm xs x(randperm(N)); % 打乱顺序破坏所有时序结构 Rs recurrenceMatrix(xs, m, tau, eps); baseOut rqaMetrics(Rs, minLen); baseline(b) baseOut.DET; end pval mean(baseline out.DET); fprintf(原始 DET %.3f, 替代数据 p %.3f\n, out.DET, pval);需要注意的是randperm 打乱会把高频噪声和低频趋势全部破坏如果原始序列本身存在强烈自相关哪怕真的没有确定性结构打乱后的 DET 也会很低这时 p 值会显得很显著。这种情况下的检验偏乐观。更严格的做法是用相位随机化对信号做 FFT把相位随机化后做 IFFT保留幅度谱但不保留相位结构。对于非线性动力学分析我给的建议是两种都跑一遍结果一致再下结论。置换检验的另一个用途是参数敏感性分析。把 m、τ、ε 各在合理范围内扫描观察 DET 的排序是否稳定。如果排序在参数微调时反复反转说明该指标对这个数据没有鲁棒性不适合做结论依据。这个检查我几乎每次都会做因为调参本身是 RQA 的最大主观来源。5.3 报告 RQA 结果时这些参数必须有名字我在看别人复现 RQA 时最常见的问题是只给一张递归图不给任何参数。递归图看起来再漂亮没有参数就等于没有定义别人无法判断你的阈值是宽松还是严格。下面是你在任何报告、笔记、实验记录里都该写清楚的参数项。参数项示例值为什么必须写数据预处理z-score 标准化不写则 ε 无法解释嵌入维 m3重构维度决定递归图分辨率延迟 τ8AMI 第一极小值延迟影响临近状态是否相关阈值 ε0.25表示“距离小于 0.25 算递归”距离范数欧氏距离换 L1 或 L∞ 指标全变最小线长 minLen2DET/LAM 口径依赖它主对角线处理已剔除影响线段长度统计有效向量数 M976边界效应和统计稳定性参考汇报时把所有参数放在一个 struct 里导出比贴在论文脚注里更可靠。我常用的做法是results.params struct(m,m,tau,tau,eps,eps,minLen,minLen,norm,euclidean,preprocess,zscore)再把 RR、DET 等指标和 params 一起存成 .mat 文件。这样即使是三个月后回来也能知道当初那个 0.83 的 DET 是怎么算出来的。6. 把 RQA 变成变化点检测器滑窗分析的一个小技巧RQA 指标最实用的工程场景不是“检验混沌”而是检测离散时间序列里的状态变化。方法很简单把序列切成固定长度的窗口在每个窗口内计算 DET 或 RR看指标沿时间轴的走势。窗口滑动到异常区域时递归结构通常会剧烈变化形成明显的低谷或尖峰。winLen 300; step 30; idxStart 1:step:length(x) - winLen 1; detSeries zeros(length(idxStart), 1); for k 1:length(idxStart) seg x(idxStart(k):idxStart(k) winLen - 1); Rk recurrenceMatrix(seg, m, tau, eps); outK rqaMetrics(Rk, minLen); detSeries(k) outK.DET; end plot(idxStart, detSeries); xlabel(窗口起始位置); ylabel(DET);这个滑窗做法的关键是窗口长度。窗口太短比如只有 80 个点重构向量太少DET 方差极大你会得到一堆假变化点窗口太长变化被平均掉检测滞后严重。我在一个模拟项目X里踩过坑窗口 200 时 DET 波动超过 0.2看起来全是变化点换成窗口 500 并做 3 点平滑才真正定位到信号突变位置。所以建议起始窗口取序列长度的 10%~20%步长不超过窗口的一半至少对比两档窗口确认变化点不是窗口参数制造出来的。使用滑窗 RQA 时阈值 ε 应该固定不能每个窗口重新用百分比阈值。否则每个窗口的 RR 都被拉齐可比较性就没了DET 的变化也会被掩盖。固定阈值后不同窗口之间的递归密度相对关系才真实。如果序列整体均值漂移明显可以在每个窗口内做 z-score 但保持 ε 固定这样兼顾归一化和可比性。那之后我养成了一个习惯任何时间序列 RQA 结果都不直接给单点值而是给出滑窗曲线和替代数据检验结果。单点的 DET 像一个快照换个窗口可能就变了滑窗曲线才真正反映结构随时间的演化也更适合跟工程故障时刻对应。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/15 17:07:43 · 阅读 1.2k RQA递归量化分析：离散时间序列的MATLAB实践与避坑指南
简介面向时间序列分析与信号处理研究者的MATLAB递归量化分析RQA工具包用于对离散时间序列构建递归图并提取复发率、最长对角线、平均线段长度和熵等量化指标适合需要考察系统非线性动态、周期性或状态切换的科研与工程场景。压缩包共2个文件包含一个可直接运行的.m脚本与一份PDF说明文档整体仅1.3MB轻量易用。已有160人学习下载。通过示例代码可快速了解数据预处理、相似性阈值设定、递归图构建到特征量化的完整流程并参考说明文档解读各指标含义从而迁移到心电图、脑电、金融波动或机械振动等实际数据上。对于具备基础MATLAB操作能力、希望掌握非线性时间序列分析方法的从业者这份代码包能提供清晰可复现的入门参考与二次开发起点。1. RQA 递归量化分析不是玄学离散时间序列能挖出多少隐藏结构某次数据分析任务里一条来自振动台架的传感器信号频谱上几乎看不出故障前后的差异功率谱密度曲线看起来都很平坦我尝试用递归量化分析RQA把离散时间序列重构成高维相空间计算递归图中对角线结构的比例故障前后 DET 指标差了近一倍问题才一下子暴露。这类代码包在 MATLAB 生态里流传很广标题里的“RQA对离散时间序列进行递归图分析.zip”就是其中一种常见实现。RQA 的核心不是画一张花哨的递归图而是把“系统轨迹有没有回到相似状态”这件事量化为一串指标适合做故障诊断、状态识别、混沌特征判断。本文按“先懂原理、再抄代码、后避坑”的顺序带你把它真正落地而不是停留在跑通 demo 的阶段。2. 递归图与 RQA 指标先搞懂 D 维空间里的“点与线”2.1 从离散序列到相空间重构嵌入维 m 和延迟 τ 怎么定离散时间序列本质上是一串标量观测值比如每隔 1 毫秒采一次电压、每个交易日记一个价格。单独看这些标量只能看到幅度变化看不出系统内部有多少自由度。RQA 的第一步就是把这串标量放到相空间里去做法是构造延迟向量X(i) [x(i), x(iτ), x(i2τ), …, x(i(m-1)τ)]其中 m 是嵌入维τ 是延迟步数i 取值范围是 1 到 N-(m-1)τ。重构后每个时刻得到一个 m 维点这些点在相空间里连成一条轨迹递归图记录的就是这条轨迹什么时候回到某个点的邻域。这个思想来自时间延迟嵌入定理实际操作中没有人会去严格验证维数而是靠指标选 m 和 τ。我一般先算平均互信息AMI随 τ 的变化取第一个局部极小值作为延迟再用假近邻FNN或 Cao 方法定 m。MATLAB 里没有内置的 AMI 函数但实现不复杂把序列分成若干等宽区间统计联合概率和边际概率计算互信息。还有更省事的替代看自相关函数取自相关第一次衰减到初始值 1/e 附近的延迟。这个做法对过采样的信号尤其有效因为过密采样会让相邻样本高度相关τ1 时递归图全是“厚对角线”好像确定性很强其实是采样率制造的假象。血泪经验是不要为了“包含更多历史信息”而把 m 设得很大。有一段 N200 的短序列我一开始用 m10重构后向量只剩 110 个距离矩阵稀疏到没法看DET 和 LAM 几乎随机跳动。后来把 m 压回 3指标才稳定。离散时间序列长度有限时m 越大有效向量越少不确定性越高我自己在工程里常用范围是 m2 到 7短序列建议从 m2、τ1 起步先看递归图形态再微调。Cao 方法的好处是它不依赖阈值比 FNN 稳定适合批处理。2.2 RQA 核心指标一览RR、DET、L、LAM、ENTR 到底在说什么有了递归矩阵 RR(i,j)1 表示第 i 个重构向量和第 j 个重构向量之间的距离小于阈值 ε。矩阵里全是点和空白人眼能看出纹理但机器要的是数字。RQA 指标就是把纹理压缩成一组统计量下面这几个是我最常用的。指标全称/计算对象高值说明RR递归点密度R 中 1 的占比状态重现的总频率过高多半是阈值太大DET大于等于 lmin 的对角线段上的点占全部递归点的比例轨迹遵循相似路径的比例代表确定性结构L所有对角线段长度的均值系统能保持相似演化路径的平均长度Lmax最长对角线段长度存在周期或很稳定的结构但可能是单根离群段LAM大于等于 vmin 的垂直线段上的点占全部递归点的比例轨迹短暂“卡”在某个区域的比例反映间歇与层流TT垂直线段平均长度系统停留在同一状态的平均时间ENTR对角线段长度分布的 Shannon 熵递归结构越复杂熵越高计算 DET 最常见的公式是DET Σ_{l≥lmin} l·P(l) / Σ_{l≥1} l·P(l)其中 P(l) 是长度为 l 的对角线段出现次数。LAM 的定义类似只不过把对角线换成垂直线。从公式能看出lmin 和 vmin 直接决定了指标口径通常在 2 到 5 之间取整数。对量化很糙的离散序列比如只取 0 和 1 的符号序列长度 1 的线段会特别多DET 会被压得很低这时候不要立刻怀疑代码先确认最小线长是否合理。需要特别提醒RQA 指标对离散序列长度非常敏感。N 在 1000 以上时 RR、DET 的估计才比较稳N 在 200 到 500 时单次计算结果波动大最好用滑动窗口或多次重复采样取平均。不要拿着 100 个点的序列去断言“系统是混沌的”那是在赌运气。2.3 阈值 ε 的选择固定阈值、百分比阈值与最大距离阈值RQA 里最容易被说成“玄学”的就是阈值 ε。“为什么你取 0.2不取 0.15”这个问题我也被问过很多次。ε 本质上是在定义“多近算相似”它必须与距离的量纲和数据动态范围一致。三种常见选法各有利弊。固定绝对阈值实现最简单但一遇到量纲不同的数据就翻车。序列 A 振幅是 0.01 量级序列 B 是 100 量级用同一个 ε 根本没有可比性。所以固定阈值之前我几乎总是先把序列标准化为 z-scorex (x - mean(x)) / std(x)。标准化之后ε 取 0.1 到 0.3 通常能给出合理密度但也要看数据形状。百分比阈值是按距离分布取分位数。比如想要递归点密度大约 15%就先算所有重构向量两两之间的距离取第 85 百分位距离作为 ε。这种做法在不同数据集之间可比性最好因为每个数据集的 RR 都被拉到同一水平。缺点是它掩盖了“原本状态重现就不频繁”的差异因此用于比较系统复杂度时没问题用于异常检测时要谨慎。最大距离比例阈值是把 ε 定义为 α 乘以距离矩阵最大值α 常见 0.05 到 0.3。这个方案不受量纲影响但最大值容易被个别离群距离带偏。我的习惯是先把距离矩阵拍平画直方图看低尾位置再综合选择。小技巧是用下面这段快速找出分位数阈值d D(:); d sort(d); idx round(0.15 * numel(d)); % 期望约15%的递归点 eps d(max(idx, 1));这里的 d 是距离矩阵展开后的向量D 来自第 3 章的递归矩阵函数。用分位数阈值有一个附带好处递归点密度不再依赖数据幅值报告结果时也更容易复现。无论选哪种请一定把 ε 和选法写进文档否则三个月后你连自己都看不懂自己算的是什么。3. MATLAB 代码落地把 RQA 跑通的最小实现与参数设置标题里的代码包解压后我认为没必要逐行去读里面的每个画图脚本先把两个核心函数理清楚一个算递归矩阵一个算 RQA 指标。我习惯把这两件事拆开因为调参时只需要换 ε、m、τ不需要碰指标计算。下面这套骨架是我在类似任务里反复用的你可以直接对照着重构自己的版本。3.1 用 MATLAB 从零写递归矩阵最小函数与内存说明function R recurrenceMatrix(x, m, tau, eps) % x: 输入离散时间序列一维列向量 % m: 嵌入维数正整数 % tau: 延迟步数正整数 % eps: 距离阈值建议在标准化后的尺度下选择 x x(:); x (x - mean(x)) / std(x); % 先做 z-score 标准化 N length(x); M N - (m - 1) * tau; % 重构向量个数 if M 2 error(数据太短重构后少于2个向量); end % 相空间重构X 的每一行是一个 m 维向量 X zeros(M, m); for i 1:M X(i, :) x(i:tau:i (m - 1) * tau).; end % 两两距离的平方利用 |a-b|^2 |a|^2 |b|^2 - 2ab X2 sum(X.^2, 2); D2 X2 - 2 * (X * X.) X2.; D2(D2 0) 0; % 消除浮点运算引入的微小负值 D sqrt(D2); R D eps; % 距离小于阈值记为递归点 R(1:M1:end) false; % 去掉主对角线自匹配 end这段代码的思路是先用矩阵运算一次性算出所有状态向量的两两欧氏距离再按阈值比较。标准化这一步不是可选的除非你明确知道数据量纲一致且要做固定阈值比较。X 的行数是 M列数是 mX * X. 是 M×M 的矩阵乘法当 M 很大时会占内存M2000 时距离矩阵约 32 MBM10000 时直接到 800 MB。所以长序列请务必分段滑动或者把距离计算改成循环逐块计算不要在一个函数里硬撑。关于阈值参数 eps注意它现在是标准化后 z-score 尺度下的绝对距离。如果前面做了标准化eps0.2 表示两个状态向量的欧氏距离小于 0.2 倍标准差。如果你更想用百分比阈值可以在调用这个函数之前算好 eps再把 eps 传进来。参数 m 和 tau 的选择见 2.1不要再让主程序里出现魔法数字。3.2 计算 RQA 指标从递归矩阵到 RR、DET、LAM、熵function out rqaMetrics(R, minLen) % R: 递归矩阵逻辑矩阵或0/1矩阵 % minLen: 最小对角/垂直线长度默认2 if nargin 2, minLen 2; end N size(R, 1); R double(R); out.RR sum(R(:)) / numel(R); % 递归点密度 % ---- 对角线段长度统计 ---- diagCounts zeros(N, 1); for d -(N-1):(N-1) seg diag(R, d) 1; L runLengths(seg); for k 1:length(L) diagCounts(L(k)) diagCounts(L(k)) 1; end end totalDiagPoints sum(diagCounts .* (1:N).); out.DET sum(diagCounts(minLen:N) .* (minLen:N).) / totalDiagPoints; out.L totalDiagPoints / max(sum(diagCounts), eps); idx find(diagCounts 0); if isempty(idx) out.Lmax 0; else out.Lmax idx(end); end p diagCounts / max(sum(diagCounts), eps); p p(p 0); out.ENTR -sum(p .* log2(p)); % 对角长度分布的信息熵 % ---- 垂直线段长度统计 ---- vertCounts zeros(N, 1); for j 1:N L runLengths(R(:, j) 1); for k 1:length(L) vertCounts(L(k)) vertCounts(L(k)) 1; end end totalVertPoints sum(vertCounts .* (1:N).); out.LAM sum(vertCounts(minLen:N) .* (minLen:N).) / totalVertPoints; out.TT totalVertPoints / max(sum(vertCounts), eps); end function L runLengths(seg) % 统计一个逻辑向量中连续1段的长度 seg seg(:); b [false; seg; false]; startIdx find(diff(b) 1); endIdx find(diff(b) -1); L endIdx - startIdx; end这个函数把递归矩阵拆成对角线和垂直线两大类结构。ripple 在于 runLengths 统计的是每个最长连续段的长度例如某条对角线上有连续 5 个点就计一条长度为 5 的线段不会重复拆成长度 1 到 5 的多段。因此 DET 的分母 totalDiagPoints 等于所有落在线段上的递归点数分子只保留长度不小于 minLen 的段上的点数两者比值不会大于 1。使用时注意 minLen 的语义它同时作用于对角线和垂直线。标准 RQA 工具里lmin 和 vmin 可以分别设置我这里为了简洁共用一个参数如果数据里孤立点特别多建议把 minLen 设成 3 再看 DET否则长度 1 的孤立点会把 DET 拉低。输出 out 是一个 struct包含 RR、DET、L、Lmax、ENTR、LAM、TT 七个字段够大多数场景使用。3.3 参数组合表不同场景下的推荐起点没有任何一组参数能通吃所有离散时间序列但下面这张表可以作为调参起点。所有场景都默认先做 z-score 标准化。数据特征mτεminLen备注长连续信号N10002~4AMI 第一极小值0.1~0.32先看递归图纹理再微调强噪声背景3~5偏小避免过采样百分比 20%3提高 DET 抗噪能力短序列N≈2002~31~3百分比 15%~30%2不要追求高精度指标明显周期信号2~4不大于主周期一半相对宽松2DET 会非常高正常现象0/1 符号序列1~210等值递归2欧氏距离不适合可用 Hamming 距离自测代码我建议用 Logistic 映射生成一段混沌序列验证整个流程是否能得到合理结果。只要 DET 明显高于同长度白噪声的 DET就说明代码链路是通的。rng(1); r 3.8; x0 0.4; N 1000; x zeros(N, 1); x(1) x0; for i 2:N x(i) r * x(i-1) * (1 - x(i-1)); % 一维混沌映射 end m 2; tau 1; eps 0.25; minLen 2; R recurrenceMatrix(x, m, tau, eps); out rqaMetrics(R, minLen); disp(out);如果 out.DET 非常低先检查阈值是否太严递归矩阵可能几乎全空如果 DET 接近 1又被对角线丛覆盖则多半是 m 过大或阈值过大。混沌映射的 DET 通常显著高于白噪声但不会接近 1这样一组自测能帮你快速定位是重构问题还是指标计算问题。4. 离散序列 RQA 避坑指南重复值、嵌入维与阈值选错的补救办法4.1 重复值导致 RR 虚高到离谱现象传感器数据经过四舍五入后出现大量重复值比如某一档位电压在序列里连续出现几十次计算出的递归图是一块块矩形色块RR 冲到 50% 以上看起来系统“天天回到同一个状态”但这其实是量化精度造成的假象。原因绝对距离阈值下完全相等的两个点距离是 0无论 ε 多小它们都算递归点。离散时间序列如果来源是数值仿真或 ADC 采样量化台阶会制造系统性重复。解决最稳妥的办法是在标准化之前检查 x 的唯一值数量如果重复值比例超过 5%要么提高数据精度要么给信号加一个幅度很小的抖动比如随机噪声振幅取数据标准差的 1e-6。抖动不会改变动力学结构但能打破完全相等的假递归。这个方法有“后悔药”的属性我一般在最后报告里注明“已加微幅扰动”避免别人质疑。4.2 嵌入维 m 过大导致递归图变成“对角线丛”现象把 m 从 2 一直调到 6递归图越来越密DET 逼近 0.95对角线又粗又长内心觉得“终于找到强确定性了”但换一段白噪声也能得到类似结果。原因m 增大后重构向量个数 MN-(m-1)τ 快速下降。向量少了距离矩阵的基数变小只要阈值稍微给大一点任意两个向量都容易“相似”于是出现伪对角线结构。解决把 M 纳入报告。一个粗糙经验是 M 至少要大于 50且 m 不要超过 N/(3τ)。更正规的是先用 Cao 方法确定 mCao 方法会输出 E1(m) 和 E2(m) 两条曲线E2 从接近 1 变成明显大于 1 的位置就是建议的 m。很多 MATLAB 实现是把 Cao 和递归图分开算我也是这么做。4.3 阈值 ε 选择不当算出来的指标全是黑匣子现象同一段数据ε 从 0.05 调到 0.5RR 从 3% 涨到 80%DET 从 0.2 涨到 0.9。把两次结果放在同一张图里结论完全相反于是有人开始把 RQA 称为“玄学”。原因ε 本质上是判断“相似”的尺度不同尺度下递归结构不同。这不是 RQA 的缺陷而是使用者没有把阈值固定下来。解决不要只算一次我习惯做 ε 扫描设置一组候选阈值计算对应的 RR 和 DET画成曲线。曲线出现明显拐点的区域通常是一个相对稳定的 ε 段如果曲线陡峭说明数据本身对尺度极敏感这时候报告结果就要特别说明。另一个可行方案是使用百分比阈值把 RR 固定住不同数据集之间就能在统一密度下比较结构差异。无论哪种方案都必须写清楚阈值策略。4.4 白噪声和周期信号都可能有较高的 DET别急着下结论现象对一段白噪声做 RQADET 也有可能到 0.1 甚至 0.2对一段理想正弦波DET 直接接近 1。于是有人看到 DET 高就声称“有确定性”翻车概率很大。原因白噪声的随机游走也能偶相遇成长度为 2 或 3 的对角线段在递归点总数不多时占比不小正弦波本身就是周期信号轨迹无限重复DET 高是必然但不能说明它“混沌”或者“具有非线性特征”。解决加一层替代数据检验。把原始序列随机打乱生成 100 组替代序列分别计算 DET看原始 DET 是否超过 95% 分位。如果连打乱后都差不多说明统计量不敏感如果明显高于才有底气说“该序列的递归结构不是偶然”。替代数据方法不完美但它能拦住最明显的过度解读。4.5 边界效应与主对角线处理方式影响指标现象递归图四个角边缘纹理明显异常主对角线附近总是有一条特别连续的亮线把主对角线置 0 之后原本长度 3 的线段被切断成两段长度 1DET 骤降。原因重构向量来自序列前部和后部时包含的样本点数量不完整主对角线表示同一个状态与自身的距离为 0如果不去除每个点都自递归也会强化对角线结构。解决报告里写明“主对角线已剔除”或“未剔除”。剔除主对角线后最小线长 minLen 通常需要设 2 或 3否则原本跨越主对角线的连续段会被误伤。处理边界效应时我一般会丢弃首尾各 (m-1)τ 个样本对应的重构向量或者至少在结果中说明计算的 M 值。很多现成工具默认剔除主对角线但如果你用的是另一套实现务必先确认它的行为。5. 让递归图说人话可视化、置换检验与参数报告规范5.1 用 MATLAB 画递归图imagesc 与 spy 的选择递归图本身是 M×M 的 0/1 矩阵最直接的可视化是用 imagesc 把矩阵画成图像。我这里会在画图前把坐标轴方向和比例尺调好否则图像上下颠倒影响阅读。figure; imagesc(1:M, 1:M, R); colormap(gray); axis xy; axis square; xlabel(状态索引 i); ylabel(状态索引 j);M 是重构向量个数。axis xy 让 y 轴从下往上递增和矩阵显示习惯一致axis square 保证图像是正方形避免递归图纹理被拉伸变形。当 M 超过 1500 时imagesc 会明显变慢这时可以用 spy 画散点对显示速度提升很大。figure; spy(R, k, 1); axis square; title(Recurrence Plot (spy view));spy 的缺点是看不出递归密度灰度只显示点位置但它适合快速浏览大矩阵。对离散时间序列来说递归图的白色区域和黑色纹理同样重要所以不要把图缩得太小否则许多细节显示不出来。保存图片时建议输出为矢量图格式因为递归图点阵缩小时会出现严重摩尔纹。5.2 用置换检验给 DET 一个可信区间RQA 指标单独看没有绝对阈值直接说“DET0.62 表示确定性很强”是不负责的。我习惯对同一段序列做一次替代数据检验用随机打乱的时间序列作为零分布。这个过程看起来简单但能挡掉很多伪结论。nPerm 100; baseline zeros(nPerm, 1); for b 1:nPerm xs x(randperm(N)); % 打乱顺序破坏所有时序结构 Rs recurrenceMatrix(xs, m, tau, eps); baseOut rqaMetrics(Rs, minLen); baseline(b) baseOut.DET; end pval mean(baseline out.DET); fprintf(原始 DET %.3f, 替代数据 p %.3f\n, out.DET, pval);需要注意的是randperm 打乱会把高频噪声和低频趋势全部破坏如果原始序列本身存在强烈自相关哪怕真的没有确定性结构打乱后的 DET 也会很低这时 p 值会显得很显著。这种情况下的检验偏乐观。更严格的做法是用相位随机化对信号做 FFT把相位随机化后做 IFFT保留幅度谱但不保留相位结构。对于非线性动力学分析我给的建议是两种都跑一遍结果一致再下结论。置换检验的另一个用途是参数敏感性分析。把 m、τ、ε 各在合理范围内扫描观察 DET 的排序是否稳定。如果排序在参数微调时反复反转说明该指标对这个数据没有鲁棒性不适合做结论依据。这个检查我几乎每次都会做因为调参本身是 RQA 的最大主观来源。5.3 报告 RQA 结果时这些参数必须有名字我在看别人复现 RQA 时最常见的问题是只给一张递归图不给任何参数。递归图看起来再漂亮没有参数就等于没有定义别人无法判断你的阈值是宽松还是严格。下面是你在任何报告、笔记、实验记录里都该写清楚的参数项。参数项示例值为什么必须写数据预处理z-score 标准化不写则 ε 无法解释嵌入维 m3重构维度决定递归图分辨率延迟 τ8AMI 第一极小值延迟影响临近状态是否相关阈值 ε0.25表示“距离小于 0.25 算递归”距离范数欧氏距离换 L1 或 L∞ 指标全变最小线长 minLen2DET/LAM 口径依赖它主对角线处理已剔除影响线段长度统计有效向量数 M976边界效应和统计稳定性参考汇报时把所有参数放在一个 struct 里导出比贴在论文脚注里更可靠。我常用的做法是results.params struct(m,m,tau,tau,eps,eps,minLen,minLen,norm,euclidean,preprocess,zscore)再把 RR、DET 等指标和 params 一起存成 .mat 文件。这样即使是三个月后回来也能知道当初那个 0.83 的 DET 是怎么算出来的。6. 把 RQA 变成变化点检测器滑窗分析的一个小技巧RQA 指标最实用的工程场景不是“检验混沌”而是检测离散时间序列里的状态变化。方法很简单把序列切成固定长度的窗口在每个窗口内计算 DET 或 RR看指标沿时间轴的走势。窗口滑动到异常区域时递归结构通常会剧烈变化形成明显的低谷或尖峰。winLen 300; step 30; idxStart 1:step:length(x) - winLen 1; detSeries zeros(length(idxStart), 1); for k 1:length(idxStart) seg x(idxStart(k):idxStart(k) winLen - 1); Rk recurrenceMatrix(seg, m, tau, eps); outK rqaMetrics(Rk, minLen); detSeries(k) outK.DET; end plot(idxStart, detSeries); xlabel(窗口起始位置); ylabel(DET);这个滑窗做法的关键是窗口长度。窗口太短比如只有 80 个点重构向量太少DET 方差极大你会得到一堆假变化点窗口太长变化被平均掉检测滞后严重。我在一个模拟项目X里踩过坑窗口 200 时 DET 波动超过 0.2看起来全是变化点换成窗口 500 并做 3 点平滑才真正定位到信号突变位置。所以建议起始窗口取序列长度的 10%~20%步长不超过窗口的一半至少对比两档窗口确认变化点不是窗口参数制造出来的。使用滑窗 RQA 时阈值 ε 应该固定不能每个窗口重新用百分比阈值。否则每个窗口的 RR 都被拉齐可比较性就没了DET 的变化也会被掩盖。固定阈值后不同窗口之间的递归密度相对关系才真实。如果序列整体均值漂移明显可以在每个窗口内做 z-score 但保持 ε 固定这样兼顾归一化和可比性。那之后我养成了一个习惯任何时间序列 RQA 结果都不直接给单点值而是给出滑窗曲线和替代数据检验结果。单点的 DET 像一个快照换个窗口可能就变了滑窗曲线才真正反映结构随时间的演化也更适合跟工程故障时刻对应。希望帮到你。本文还有配套的精品资源点击获取