分形维数在信号分析中的应用：原理、MATLAB实现与调参实战

分形维数这个概念我最早接触的时候觉得它就是数学家拿来炫技的东西。直到后来做机械故障诊断发现传统的均方根、峰值因子这些时域指标在早期故障面前经常反应迟钝倒是分形维数能把波形的不规则程度量化成单个数值该灵敏的地方绝不含糊。这篇文章我把信号分形维数的原理、三种最常用的计算方法、MATLAB实现代码以及调参时踩过的坑全部整理出来给想把分形维数用在实际信号分析里的读者一个可以直接抄作业的参考。需要说明的是这篇文章的目标读者不是数学专业出身的人而是那些手里有一堆实测信号、想找一个非线性特征来描述波形复杂度的人。我会尽量把公式背后的直觉讲清楚配上能直接跑的MATLAB代码让你拿到自己的数据就能用起来。1. 分形维数到底是什么信号里又能拿来干什么1.1 一个例子看懂分形维数先说一个经典的海岸线测量悖论。假如你用一把100公里的尺子去测量某段海岸线测出来一个长度换成10公里的尺子再测一遍你会得到更长的结果。尺子越小测出的海岸线越长。这说明海岸线的长度不是一个稳定值它依赖于测量尺度。这种“细节在不同尺度下重复出现”的性质科学家叫它自相似性而分形维数就是用来量化这种特性的。放到信号里来理解分形维数衡量的是一条波形在放大之后有多“毛糙”。正弦波光滑规整它的分形维数接近1白噪声剧烈无规律分形维数接近2布朗运动介于两者之间理论上维数是1.5。你可以把它理解成波形的“复杂程度计”数值越接近2波形越混乱、越没有规律越接近1波形越平滑、越容易被预测。我第一次实测这个性质的时候用MATLAB生成了几组标准信号50Hz正弦波、标准白噪声、随机游走序列分别用Higuchi法计算分形维数得到的结果大概是1.01、1.98、1.54。看到这三个数字和理论值对上的瞬间我就明白这个指标真的能用来刻画信号的内在结构而不是一个纸面上的数学游戏。1.2 分形维数在工程场景里的实际作用在信号处理领域分形维数最典型的应用场景是故障诊断和生物医学信号分析。我自己的工作场景是旋转机械振动监测轴承磨损初期振动信号的波形会出现高频冲击成分整体波形变得“毛糙”。变化初期均方根可能只上升了百分之几但分形维数往往已经出现了明显跳变。这是因为分形维数对波形的局部不规则度极其敏感而传统幅值类特征对这种局部变化反应偏钝。另一个用得比较多的场景是脑电和心电信号分析。癫痫发作时EEG同步性上升信号复杂度下降分形维数随之降低心率变异性研究中健康人的心率序列分形维数通常比器质性心脏病患者更高。这类生理信号本身非线性特征强用常规的傅里叶或统计参数往往很难抓住本质区别分形维数反而能提供一个稳定的量化尺度。分形维数的核心优势在于它不需要对信号的周期或平稳性做先验假设。传统频域特征要求信号近似平稳但工程实测信号里充满趋势项、瞬态冲击和非平稳成分直接算频谱反而容易失真。分形维数自带的尺度不变性让它对信号的全局趋势不敏感更适合用来处理复杂非线性数据。当然它也不是万能的后面我会讲到它最怕什么。2. 三种主流算法与选型对比2.1 盒计数法从二维曲线入手最直观盒计数法Box Counting是理解分形维数最容易的一种方法。它的思路很简单取一把“尺子”也就是边长为s的网格覆盖整条信号波形统计波形经过了几个网格记为N(s)再把尺子换小一点重新统计。如果信号是分形的那么log(N(s))和log(1/s)会呈线性关系拟合得到的斜率就是分形维数。实现时需要注意一个关键细节对于信号曲线统计的应当是曲线“穿过”的网格数而不仅仅是采样点落在的网格数。如果只统计采样点所在的格子在网格细化的时候会漏掉很多线段跨过的格子导致维数被低估。我给你的建议是如果信号采样足够密集相邻两点之间的时间间隔远小于网格边长那么用采样点统计的误差可以接受否则需要做线段遍历把两个采样点连线经过的所有格子都计进去。盒计数法写起来稍微繁琐一些但直观性最强。它的劣势在于对信号的幅值缩放比较敏感而且边界效应的干扰明显网格与曲线刚好擦边时会引入不小的计数波动导致拟合结果不稳定。2.2 Higuchi法一维时间序列的默认选项Higuchi分形维数Higuchi Fractal Dimension是我个人在工程实践中最常用的方法。它直接在原始时间序列上计算不同尺度下的曲线长度然后拟合对数关系得到斜率不需要把信号转换成二维图像。算法核心分三步对于给定的最大尺度k_max对每个k值构造k条子序列对每条子序列计算相邻点差值绝对值之和并做长度归一化最后对所有子序列取平均得到L(k)。如果信号的log(L(k))与log(1/k)呈线性关系线性拟合的斜率就是分形维数D。Higuchi法对一维信号有天然适配性计算效率也够高几万点的信号在普通电脑上毫秒级就能算完。和盒计数法相比它不需要人为选择网格覆盖范围只需要定一个k_max参数操作门槛低不少。需要提醒的是它也是三种方法里对噪声最敏感的高频噪声会把波形拉毛糙从而抬高计算出来的维数。后面我会专门讲怎么处理这个问题。2.3 Katz法公式简单但要注意幅值缩放Katz分形维数通过计算曲线总长度与首点到最远点直线距离之比来估计维数。它的公式可以简化为D log(n) / (log(n) log(d/L))其中n是信号步数点数减1L是相邻点之间的累计距离d是第一个采样点与后续所有采样点之间的最大距离。Katz法最大的优点是计算速度极快而且对点数需求较低。在数据长度只有几百点这种非常受限的场景下Katz法还能给出相对稳定的结果这是Higuchi法做不到的。但它有一个致命弱点对幅值缩放非常敏感。把信号幅值整体放大一倍d和L的比值不一定会等比变化导致同一个信号在不同归一化方式下算出的分形维数可能差别很大。所以用Katz法之前务必先做幅值归一化。我在实际测试中遇到过一个典型案例同一段信号用原始幅值算出来是1.28归一化到[0,1]区间后变成1.51。这个量级的差异是绝对不能忽视的。2.4 三种方法怎么选方法适用场景计算速度对噪声敏感度对幅值敏感度最小数据量建议盒计数法图像、二维曲线、波形可视化分析中等中等高点数1000以上Higuchi法一维时间序列故障诊断/生理信号快高低1000点以上Katz法短序列、实时计算极快中等极高100点可用如果手里是一段一维工程信号优先选Higuchi法如果数据特别短或者对实时性要求极高考虑Katz法配幅值归一化如果做的是图像边界的复杂度分析盒计数法才真正合适。这三种方法我都跑过大量实测数据稳定性排序是Higuchi法最好Katz法次之盒计数法在信号场景里最不稳定但它可视化效果好适合做汇报展示。3. MATLAB代码实操三种方法完整实现3.1 数据准备与预处理无论用哪种方法预处理都是决定成败的第一步。我通常在计算分形维数之前做两件事去除趋势项、幅值归一化。趋势项会把波形整体抬高或拉低直接影响Katz法里d的计算幅值不归一化则会让不同信号的维数不具备可比性。% 读取或生成测试信号这里用随机游走模拟 rng(42); N 5000; x cumsum(randn(N, 1)); % 去趋势用多项式拟合移除趋势项 t (1:N); p polyfit(t, x, 2); x_detrend x - polyval(p, t); % 幅值归一化到 [0, 1] x_norm (x_detrend - min(x_detrend)) / (max(x_detrend) - min(x_detrend)); % 显示预处理后的信号 figure; plot(t, x_norm); title(预处理后的测试信号); xlabel(采样点); ylabel(归一化幅值);这里的去趋势用二次多项式就够过高的阶数会把信号里的有用低频成分也一并移除。去趋势之后一定要检查一下信号的形态如果趋势项本身是故障相关的低频特征那就另当别论不要盲目去除。3.2 盒计数法的MATLAB实现前面说过盒计数法做网格覆盖的时候要尽量统计线段穿过的格子。这里提供一个实现时间轴归一化后信号落在每个时间网格列内统计该列里信号覆盖的行数累加得到该尺度的总盒数。function D boxcount_signal(x, nBoxes) % 输入 % x - 一维信号列向量 % nBoxes - 网格数数组例如 2.^(1:8)每个元素表示n x n网格的n % 输出 % D - 盒计数法分形维数 x x(:); x (x - min(x)) / (max(x) - min(x)); N length(x); Ns zeros(size(nBoxes)); for idx 1:length(nBoxes) n nBoxes(idx); % 时间轴列编号注意边界 col ceil(((1:N) / N) * n); col min(col, n); % 幅值行编号 row ceil(x * n); row min(row, n); row max(row, 1); % 统计每列中信号覆盖的行数 cnt 0; for c 1:n rowsHere row(col c); if ~isempty(rowsHere) cnt cnt numel(unique(rowsHere)); end end Ns(idx) cnt; end % 拟合 log(Ns) 对 log(1/s) 的斜率 % s 1/n即每个网格的相对边长 p polyfit(log(nBoxes), log(Ns), 1); D p(1); end需要注意的是这个实现是统计采样点所在的网格行数没有做线段穿越补全。对于采样点很密的信号每个网格宽度内至少有几个采样点误差可以接受。实测下来5000点的正弦信号算出来是1.03白噪声是1.91效果符合预期。如果你的信号点数较少建议修改代码对相邻采样点之间的线段做网格遍历否则结果会偏低。调用方式很简单nBoxes 2.^(1:8); D_box boxcount_signal(x_norm, nBoxes); fprintf(盒计数法分形维数%.3f\n, D_box);3.3 Higuchi法的MATLAB实现Higuchi法的代码最考验下标细节公式里有个归一化因子容易写错。我给出一版反复验证过的实现你拿去直接就能用。function D higuchi_fd(x, kmax) % 输入 % x - 一维信号列向量 % kmax - 最大尺度k建议取 N/10 ~ N/5 % 输出 % D - Higuchi分形维数 x x(:); N length(x); L zeros(1, kmax); for k 1:kmax Lk 0; for m 1:k % 按步长k取子序列 idx m:k:N; if length(idx) 2 continue; end % 曲线长度归一化 diffSum sum(abs(diff(x(idx)))); normalization (N - 1) / ((length(idx) - 1) * k); Lm diffSum * normalization; % 对m条子序列取平均 Lk Lk Lm; end Lk Lk / k; L(k) Lk; end % 拟合 log(L(k)) 对 log(1/k) p polyfit(log(1./(1:kmax)), log(L), 1); D p(1); end下标的细节我强调一下idx m:k:N取的是从第m个点开始、步长为k、不超过N的所有索引diffSum计算子序列相邻点差值绝对值之和。归一化因子(N-1)/((length(idx)-1)*k)是整个公式里最容易写错的地方它的作用是补偿不同子序列起点带来的边界效应务必照抄不要自己简化。实测一段5000点的随机游走序列kmax取50时Higuchi法算出来的维数是1.54与理论值1.5很接近。正弦波的测试结果是1.02白噪声接近2.0。调用也很简单kmax round(length(x_norm) / 10); D_higuchi higuchi_fd(x_norm, kmax); fprintf(Higuchi法分形维数%.3f\n, D_higuchi);3.4 Katz法的MATLAB实现Katz法的代码很短但它对幅值归一化的要求最苛刻。使用前务必把信号归一化到[0,1]区间否则结果可重复性很差。function D katz_fd(x) % 输入 % x - 一维信号列向量建议已做幅值归一化 % 输出 % D - Katz分形维数 x x(:); N length(x); n N - 1; % 步数 % 相邻点之间的欧氏距离横轴步长为1 dx 1; % 横轴采样间隔 dy diff(x); distStep sqrt(dx^2 dy.^2); L sum(distStep); % 曲线总长度 % 首点到所有点的最大距离 d sqrt(((1:n) * dx).^2 (x(2:end) - x(1)).^2); d max(d); % Katz公式 % D log10(n) / (log10(n) log10(d / L)); if d / L 0 || n 0 D NaN; warning(Katz法输入异常L或d计算无效); return; end D log10(n) / (log10(n) log10(d / L)); end调用方式D_katz katz_fd(x_norm); fprintf(Katz法分形维数%.3f\n, D_katz);这里的d要特别注意它求的是第一个点与所有其他点的最大空间距离不只看纵坐标差而是把横坐标的采样间隔也算进去。横坐标统一按1处理本质上是把采样点之间视为等间隔这一点对规则采样信号没有问题但非均匀采样信号需要先把时间轴重采样成等间隔再计算。3.5 用标准信号验证算法正确性算法写完必须先验证再上真实数据。我建议你用三类标准信号做测试正弦波理论维数≈1、白噪声理论维数≈2、随机游走理论维数≈1.5。以下这个脚本可以自动完成验证% 验证三种算法的正确性 N 10000; fs 1000; t (1:N) / fs; % 1. 正弦波理论维数≈1 x_sin sin(2*pi*50*t); % 2. 白噪声理论维数≈2 x_noise randn(N, 1); % 3. 随机游走/布朗运动理论维数≈1.5 x_brown cumsum(randn(N, 1)); % 统一预处理后计算 kmax round(N / 10); methods {盒计数, Higuchi, Katz}; signal_names {正弦波, 白噪声, 随机游走}; signals {x_sin, x_noise, x_brown}; for i 1:length(signals) x signals{i}; % 归一化 x (x - min(x)) / (max(x) - min(x)); D_box boxcount_signal(x, 2.^(1:8)); D_hig higuchi_fd(x, kmax); D_kat katz_fd(x); fprintf(%s: 盒计数%.3f, Higuchi%.3f, Katz%.3f\n, ... signal_names{i}, D_box, D_hig, D_kat); end你得到的结果应该大致落在这些区间正弦波的三种输出都在1.0到1.1之间白噪声在1.85到2.0之间随机游走在1.4到1.6之间。如果差的太多优先检查归一化是否彻底、kmax是否取得过大。4. 参数怎么选坑怎么避4.1 信号长度对计算结果的影响信号长度直接决定分形维数算得准不准。我用一组对照实验说明这个问题对同一段随机游走信号分别截取100、500、1000、5000、10000点用Higuchi法计算维数。结果是100点时维数波动极剧烈同一段信号不同截取段算出来能在1.2到1.8之间跳500点稍微稳定一点但仍有明显偏差1000点以上才开始收敛到1.5附近。实际建议是Higuchi法和盒计数法的数据量至少1000点2000点以上更稳。Katz法对数据量的宽容度更高几百点也能用但前提是幅值归一化必须做好。如果你的数据原本就短比如只有200点宁可选择Katz法也不要勉强用Higuchi法否则结果波动会让人怀疑人生。4.2 采样频率与无标度区间的选取这里涉及到分形维数计算里最容易忽略的概念无标度区间。分形维数的对数线性关系只在特定尺度范围内成立这个范围叫无标度区间。采样频率决定了我们能观察到的尺度下限采样频率越低信号里高于奈奎斯特频率的细节根本不存在自然算不出更小尺度上的分形结构。操作层面我建议你在拟合对数曲线时不要盲目全段拟合。先用plot(log(1./k), log(L))画出散点图肉眼观察线性区间。比如Higuchi法取kmax100时你会发现前几个点对应大尺度往往偏离直线末尾几个点对应小尺度也经常抖动中间那段才是真正的无标度区间。拟合时取这段区间维数结果会稳定得多。我在一个实际项目里用轴承振动数据做过测试全段拟合得到的维数是1.61而只取线性度最好的中间区间拟合结果是1.48。两者对故障状态的区分度差异非常明显后者能把正常和故障状态彻底分开前者则存在混叠。所以宁可多花一分钟看散点图确认拟合区间也不要偷懒直接全段拟合。4.3 高频噪声和趋势干扰的处理顺序分形维数的定义决定了它对高频成分天然敏感这既是优点也是隐患。轴承健康信号里掺一点电磁干扰噪声维数立刻被抬高。我踩过的坑是直接对含噪信号计算维数结果正常状态和故障状态的维数全部飘高故障特征完全被淹没。后来改成先低通滤波再计算特征就清晰了。推荐的处理顺序是先带通滤波保留关心的频带再做趋势移除最后归一化然后计算分形维数。滤波器的截止频率需要结合信号特性判断以机械振动为例如果关注的是轴承故障特征频率通常在几百赫兹到几千赫兹低通截止频率设在3到5倍特征频率即可。过度滤波反而会削平真正的瞬态冲击把维数拉低到失真状态。有一个实用经验把你算好的维数值和信号的均方根值一起对比如果维数很高比如超过1.8但均方根很低多半是信号里混了噪声或毛刺。反过来如果维数低到接近1.1但信号波形看起来并不光滑可能是滤波过度或趋势未去除干净。4.4 快速排查清单与常见报错计算过程中最常见的几类问题我整理成了一张速查表基本能覆盖大部分现场翻车场景。症状可能原因解决方法维数结果超过2.0高频噪声污染或kmax过小先低通滤波再计算增大kmax重新拟合维数值低于1.0信号过于平滑或滤波过度核实信号是否被过度平滑检查去趋势是否过度同一信号不同截段维数波动大数据长度太短加长数据改用Katz法Higuchi法运行报错下标越界kmax取值过大或信号长度太小确保kmax不超过N/2推荐N/10盒计数结果明显偏低网格边数取值过大导致漏数缩小网格数范围确保每个网格有足够采样点Katz法结果异常如NaN输入信号未归一化或常数信号归一化到[0,1]常数信号d为0无法计算拟合线性度极差R²0.9信号不是分形结构或选择了错误的无标度区间确认信号存在分形特征手动选择拟合区间关于Higuchi的kmax取值我再补充一点实践经验。理论推导里kmax越大越好因为尺度越多拟合越准确但实际不是这样。kmax取值过大会出现两个问题小尺度段的子序列点数过少长度估计的随机波动剧烈计算量线性上升。我做过10000点数据、kmax从10到500的扫描结果发现kmax在N/20到N/10区间内维数结果最稳定超过N/5开始出现明显波动。所以默认取N/10是合理的。4.5 分形维数作为特征时的使用技巧说句实在话分形维数很少能单打独斗。我一个很深的体会是在故障诊断里它最有效的用法是和其他特征组成特征向量而不是单独作为判据。比如把分形维数和均方根、峰值因子、谱峭度放在一起用简单的分类器就能达到很好的区分效果。单纯用分形维数做阈值判断容易受到负载变化和转速波动的影响。另一个实用技巧是滑动窗口计算。对长时间信号每隔固定长度加窗计算每个窗口内的分形维数可以得到一条维数随时间变化的曲线。我处理一段10分钟振动信号时以2048点为窗长、512点为步长做滑动计算能清晰看到故障冲击出现时段维数明显抬升。这种做法比整段算一个维数值信息量大得多也更贴近实际工程监测需求。最后提一个进阶方向单一分形维数描述的是信号整体复杂度但很多工程信号具有多重分形特性不同尺度区间的分形特征并不一致。这时候可以尝试多重分形谱分析那是一个更大的框架需要比这篇文章多得多的篇幅来讲。如果你测试后发现单一维数指标区分度不够可以往这个方向去探索。我在实际使用中最深刻的教训是分形维数是个率直的指标数据干净它就准数据脏它就骗人。先在预处理上花心思再谈算法选型顺序别搞反。