简介这份资源聚焦Matlab环境下的多重分形谱算法与盒子维数计算面向分形几何、图像处理、信号分析及金融建模等方向的学习者与研究者帮助其用数值手段定量刻画复杂系统的自相似性与尺度不变性。压缩包共2个文件均为.m脚本整体约3KB分别对应多重分形谱算法与盒子分形维数的实现可直接在Matlab中运行调试。已有1343人学习下载说明其在分形分析入门与实践中具有一定参考价值。内容围绕数据分箱、盒计数、H(q)与f(α)特征谱构造以及log-log图拟合求盒子维数等核心环节展开读者可据此理解多重分形谱的完整计算流程掌握盒子维数的参数设定与斜率拟合思路并在此基础上迁移到图像、时间序列等实际数据的分析中适合作为算法复现与课程实验的参考脚本。1. 多重分形谱与盒子维数为什么一条曲线能暴露数据的“不均匀程度”很多人第一次接触分形维数会下意识觉得它就是个“小数版本的维度”——直线是 1平面是 2那 1.7 大概就是某种介于两者之间的形状。这个直觉不算错但它只描述了分形最粗的那一层。真正让分形在信号分析、图像处理、材料表征里站住脚的是另一件事同样一条曲线它的“不均匀程度”是可以被量化的而多重分形谱就是干这个的。举个反直觉的例子。两条长度完全相同的时序曲线一条处处平滑一条在某些局部剧烈抖动、其余地方几乎不动。它们的盒子维数可能非常接近甚至一样。但多重分形谱一算差别立刻暴露剧烈抖动的那条谱宽明显更大。这说明盒子维数只告诉你“整体有多占空间”而多重分形谱告诉你“不同区域各自有多占空间”。前者是一个数后者是一条曲线。这就是为什么这套算法值得搞懂。盒子维数是入口多重分形谱是真正的分析工具。很多做图像纹理、故障诊断、金融波动分析的人卡就卡在只会算一个维数却解释不了“为什么两类样本维数差不多但性质完全不同”。多重分形谱解决的正是这个解释力问题。接下来我会从盒子维数的计数机制讲起再推到多重分形谱的配分函数和 Legendre 变换最后落到 MATLAB 里怎么把这条谱稳定地算出来。2. 盒子维数从“数盒子”到覆盖计数机制到底怎么跑2.1 盒子计数法的核心动作换尺度、数非空格子盒子维数的定义看起来简单用边长为 ε 的格子去覆盖目标集合数出有多少个格子至少包含集合的一个点记作 N(ε)。然后不断缩小 ε看 N(ε) 怎么增长。如果满足 N(ε) ∝ ε^(-D)那 D 就是盒子维数。关键在于“怎么数”。对一条离散曲线你不能真的拿几何盒子去套而是把数据映射到网格上。常见的实现思路是先把曲线归一化到一个单位正方形然后按当前 ε 划分网格遍历所有数据点标记它落在哪个格子最后统计被标记过的格子数。这里有个容易翻车的点——曲线是连续的但数据点是离散的如果采样点太稀小尺度下会漏掉很多本该被覆盖的格子导致 N(ε) 偏小维数被低估。所以盒子维数的计算精度本质上取决于两件事数据点密度是否足够支撑到你想要的最小尺度以及尺度序列怎么选。尺度不能随便取通常取 2 的幂次或者等比例序列保证 log(1/ε) 在横轴上分布均匀拟合时才不会因为某一段点太挤而带偏斜率。2.2 双对数拟合为什么斜率就是维数把 log N(ε) 对 log(1/ε) 做线性回归斜率就是盒子维数。这一步的机制来自定义本身N(ε) ∝ ε^(-D) 两边取对数得到 log N(ε) D · log(1/ε) C。所以拟合斜率直接对应 D。但实际做的时候不是所有尺度都可用。ε 太大时格子数太少统计涨落大ε 太小时离散数据不够N(ε) 饱和不再增长。这两个区间都会让双对数曲线弯曲。合格的实现应该只取中间那段近似线性的区间做拟合或者至少把拟合残差输出出来让人判断哪些尺度被污染了。% 盒子维数核心多尺度覆盖计数 双对数拟合 function [D, logN, logInvEps] boxDim(points, epsList) % points: N×2 归一化坐标epsList: 尺度序列从大到小 logN zeros(size(epsList)); logInvEps zeros(size(epsList)); for k 1:length(epsList) eps epsList(k); % 把点映射到网格索引唯一化后统计非空格子数 idx floor(points / eps) 1; occupied unique(idx, rows); logN(k) log(size(occupied, 1)); logInvEps(k) log(1/eps); end % 最小二乘拟合斜率即盒子维数 p polyfit(logInvEps, logN, 1); D p(1); end这段代码里floor(points/eps)1是网格索引映射unique(...,rows)完成非空格子去重计数。polyfit做一次线性拟合斜率p(1)就是维数。参数epsList的选择直接决定结果稳定性我一般会取从 1/4 到 1/64 左右的一组等比尺度避开两端饱和区。如果拟合残差平方和偏大说明尺度区间没选好需要重新截取。2.3 盒子维数的局限一个数装不下不均匀性盒子维数把所有格子一视同仁只关心“有没有点”不关心“有多少点”。这就导致它对局部密度差异不敏感。一条曲线上如果有一小段特别密集、其余稀疏盒子维数几乎只反映稀疏部分的覆盖行为密集段的信息被平均掉了。这个局限不是实现问题是定义本身的天花板。要突破它就必须引入“测度”的概念——不再只数格子有没有点而是看每个格子里有多少点按密度加权。这就是多重分形谱的起点。3. 多重分形谱配分函数如何把“密度差异”拆成一条曲线3.1 从计数到加权配分函数的引入多重分形的核心动作是把“格子里的点数”变成一个概率测度。设第 i 个格子的点数为 n_i总点数为 N则 p_i n_i / N。然后构造配分函数χ(q, ε) Σ p_i^q这里的 q 是一个实数阶矩。q 取不同值相当于用不同力度放大高密度区域和低密度区域。q 很大时p_i^q 会被大 p_i 主导配分函数只反映密集区q 很小时小 p_i 的贡献被放大反映稀疏区。q 0 时p_i^0 1配分函数退化成格子数回到盒子维数的情形。这就是多重分形“多重”的含义通过扫描 q把不同密度层次的标度行为分开。如果存在标度关系 χ(q, ε) ∝ ε^τ(q)那 τ(q) 就刻画了第 q 阶的标度指数。对 τ(q) 做 Legendre 变换就得到多重分形谱 f(α)α dτ/dqf(α) q·α - τ(q)α 叫奇异指数描述局部密度奇异性f(α) 是具有该奇异指数的子集的分形维数。谱越宽说明密度分布越不均匀。3.2 配分函数的数值陷阱q 扫描与归一化理论很干净数值实现却容易翻车。第一个坑是 p_i 可能为 0。如果某个格子在当前尺度下没有点p_i 0那 q 0 时 0^q 会变成 Inf直接把配分函数炸掉。常见的处理是只对非空格子求和或者给 p_i 加一个极小量。但加极小量会改变小 q 段的行为所以我更倾向只对非空格子计算并在 q 为负时确认非空格子数足够多。第二个坑是 q 的步长和范围。q 太密计算量大q 太疏Legendre 变换出来的谱会有锯齿。一般取 q 从 -10 到 10步长 0.5 左右能覆盖大多数实际信号的谱宽。如果谱在两端还没收敛就截断f(α) 会出现假的开口。% 配分函数与 tau(q) 计算 function [tau, qList] partitionFunction(points, epsList, qList) tau zeros(size(qList)); for qi 1:length(qList) q qList(qi); logChi zeros(size(epsList)); for k 1:length(epsList) eps epsList(k); idx floor(points / eps) 1; [~, ~, ic] unique(idx, rows); counts accumarray(ic, 1); % 每个格子的点数 p counts / sum(counts); % 归一化为测度 if q 1 logChi(k) sum(p .* log(p)); % q1 时取极限形式 else logChi(k) log(sum(p.^q)); end end % 对 logChi 与 log(eps) 拟合斜率得到 tau(q) pf polyfit(log(epsList), logChi, 1); tau(qi) pf(1); end end这里accumarray统计每个格子的点数p counts/sum(counts)完成测度归一化。q 1 时sum(p.^q)恒为 1取对数后为 0无法提供标度信息所以单独用sum(p.*log(p))的极限形式。polyfit对每个 q 拟合出 τ(q)。参数qList的范围要覆盖谱的两端否则 Legendre 变换后谱形不完整。3.3 Legendre 变换从 τ(q) 到 f(α) 的最后一跳有了 τ(q) 序列下一步是数值微分求 α dτ/dq再算 f(α) q·α - τ(q)。数值微分的稳定性很关键直接差分在 q 步长较大时会放大噪声。我一般用中心差分或者先对 τ(q) 做一次平滑再求导。% Legendre 变换tau(q) - f(alpha) function [alpha, f] legendreTransform(qList, tau) alpha gradient(tau, qList); % 中心差分求 dτ/dq f qList .* alpha - tau; % 计算 f(α) endgradient默认用中心差分比手动diff更稳。得到的alpha和f就是多重分形谱的横纵坐标。谱宽 Δα max(α) - min(α) 是衡量不均匀性的关键指标Δα 越大分布越不均匀。如果算出来的 f(α) 出现负值或者明显不对称的尖刺通常是 τ(q) 拟合区间没选好或者 q 范围不够。4. 避坑与排查这套算法最容易翻车的几个地方4.1 现象双对数曲线明显弯曲拟合残差大原因通常是尺度序列跨越了饱和区。ε 太大时格子数少、统计涨落大ε 太小时离散点不够、N(ε) 不再增长。这两个区间都会让 log N 对 log(1/ε) 偏离直线。解决方式是先画出双对数散点图肉眼确认线性区间再截取中间段做拟合。不要一上来就把所有尺度丢进polyfit。我一般会输出每个尺度的 logN 值检查末尾几个点是否已经持平。4.2 现象q 为负时配分函数出现 Inf 或 NaN原因是某些格子在当前尺度下点数为 0p_i 0而 0 的负次幂是 Inf。如果代码里直接对全网格求和就会污染整个配分函数。解决方式是只对非空格子求和并在 q 0 时检查非空格子数是否足够。如果非空格子太少说明当前尺度已经低于数据支撑能力应该把该尺度从负 q 的计算中剔除或者整体缩小尺度范围。4.3 现象多重分形谱出现锯齿或负值原因通常是 τ(q) 的数值微分不稳定。q 步长太大、τ(q) 本身有噪声、或者 q 范围两端没收敛都会让gradient出来的 α 抖动进而让 f(α) 出现负值。解决方式是减小 q 步长、对 τ(q) 做轻度平滑、或者扩大 q 范围直到谱两端闭合。如果平滑后谱形仍然异常要回头检查配分函数拟合的线性度而不是继续在微分环节调参。4.4 现象两类样本盒子维数接近但多重分形谱差异明显这不是 bug是正常现象但容易被误判为算法不稳定。盒子维数只反映覆盖计数对密度差异不敏感多重分形谱通过 q 加权把密度层次拆开所以能区分。解决方式是在报告结果时同时给出盒子维数和谱宽 Δα不要只报一个数。如果业务上关心的是局部不均匀性盒子维数只能作为辅助主指标应该是 Δα 和谱形对称性。4.5 现象换一组数据后谱形完全变了怀疑算法不可复现原因往往是归一化方式不一致。不同数据的坐标范围不同如果没统一归一化到同一尺度网格划分的物理含义就变了N(ε) 和 p_i 都会漂移。解决方式是在计算前强制把数据归一化到单位正方形或单位区间并固定尺度序列。归一化方式要在文档里写清楚否则同一套代码换个人跑结果对不上。5. 设计取舍为什么这套算法值得自己实现一遍把盒子维数和多重分形谱放在一起看会发现它们其实是一套递进的设计盒子维数用“有没有点”回答“占不占空间”多重分形谱用“有多少点”回答“怎么占空间”。前者是后者的 q 0 特例后者是前者的密度加权推广。理解了这个关系再去看 MATLAB 里的实现就不会把两段代码当成孤立工具而会明白它们共享同一套网格计数和双对数拟合骨架。这套算法真正的难点不在公式而在数值细节的取舍。尺度序列选多长、q 范围取多宽、拟合区间怎么截、负 q 怎么处理每一个选择都会影响最终谱形。这些取舍没有唯一正确答案取决于数据特性和分析目的。比如做故障诊断时谱宽 Δα 的稳定性比绝对精度更重要那就可以牺牲一点分辨率换平滑做材料表征时谱形对称性可能更关键那就需要更密的 q 扫描和更严格的拟合区间控制。我自己的习惯是每次算完先画三张图双对数覆盖计数图、τ(q) 曲线、f(α) 谱。三张图对不上就说明中间某一步的参数有问题而不是直接信最后那个数。这个习惯帮我避开了很多“看起来算出来了、其实不可信”的结果。从那以后我每次拿到新数据都强制走一遍这三张图的检查流程再决定要不要调尺度或 q 范围。如果你手头正好有需要做分形分析的曲线或图像数据这套 MATLAB 实现可以直接拿来跑也可以顺着上面的机制说明改尺度策略和 q 扫描逻辑。希望帮到你。本文还有配套的精品资源点击获取