简介多重分形谱算法与盒子维数计算是分形几何中分析复杂系统自相似结构的重要工具这套基于Matlab编写的代码包面向需要量化研究多维数据特征的科研人员和工程技术人员可直接用于数值实验与教学演示。压缩包内共包含2个m文件分别实现多重分形谱的完整计算流程与盒子分形维数的自动求解压缩后大小仅为3KB极为轻量便于下载、查看与运行。脚本覆盖了从数据预处理、分箱统计、盒计数到特征谱构造、log-log图线性拟合等关键环节代码结构清晰注释友好适合具备基础Matlab知识的学习者快速上手。资源上线以来已有1341人学习或下载获得一定关注与使用反馈。通过运行这两个脚本读者不仅能够理解多重分形谱与盒子维数背后的数学原理还能掌握相关的编程实现技巧并迁移到图像纹理分析、信号奇异性检测、金融市场波动建模等实际场景中为后续研究提供有力支撑。 做信号处理和图像分析的朋友大概率都碰到过“分形维数”这个词。我在一个材料微观结构分析项目里需要量化孔隙复杂程度查文献后发现相关方法基本都靠MATLAB核心就是盒子维数和多重分形谱算法。当时网上的资料要么只给代码不讲原理要么公式堆了一大堆却和程序对不上号我折腾了蛮久才把完整流程跑通。这篇文章把两个算法的原理、MATLAB实现、参数设置和常见坑一次讲清楚适合刚装好MATLAB、想入门分形分析以及被D0、D1、f(α)这些符号绕晕的朋友。看完之后你能直接动手算出一张图像或一段信号的盒子维数和多重分形谱也能明白代码里每一步在干什么改参数时心里有底。1. 先理清概念盒子维数与多重分形谱是什么关系1.1 为什么非整数维度能描述复杂几何在经典几何里点是0维、直线是1维、平面是2维、立方体是3维这些都是整数维度。但自然界中的海岸线、裂纹、云层边界、血管网络远不止这么简单。用一把固定长度的尺子去量海岸线尺子越短测量到的长度越长这说明海岸线的真实“维度”无法用1维或2维简单描述。Mandelbrot把这种介于整数之间的维度称为分形维数用来刻画几何对象的复杂程度和自相似性。分形维数越大说明对象在有限空间里填充得越“满”细节越丰富。1.2 盒子维数工程上最实用的分形维数算法理论上有Hausdorff维数等严格定义但工程上计算不便。盒计数法box-counting method成为最常用的近似方法它又被称为盒子维数。把目标图像划分成边长为δ的网格统计至少包含一个目标点的格子数量N(δ)然后不断改变δ观察N(δ)的变化。如果对象是分形的N(δ)与1/δ在双对数坐标下会呈现线性关系拟合斜率就是盒子维数D。在项目里我会把孔隙图像二值化后直接套用盒计数法图像分辨率决定了δ的下限这直接影响到最终维数的精度。1.3 多重分形谱一个维度数解决不了的问题盒子维数描述的是整体复杂度但对很多实际对象不同区域或不同尺度下的分形行为并不一样。比如岩石孔隙中有些区域孔隙密集、连通强有些区域稀疏孤立在这种情形下单一维数会把不同局部的特征平均掉了。多重分形谱算法不是返回一个数而是返回一条曲线横坐标α称为奇异强度指数纵坐标f(α)表示奇异强度为α的所有子集的分形维数。通过这条曲线可以了解测度在不同位置的分布是否均匀、差异有多大。当对象是均匀分形时f(α)曲线会收缩成一个点这个点的纵坐标正好等于盒子维数D0。可以说盒子维数是多重分形谱的一个特例多重分形谱则是盒子维数在更复杂对象上的延伸。2. 算法原理拆解从盒计数到多重分形谱2.1 盒子维数怎么算公式、步骤、拟合细节盒子维数的数学定义是D lim(δ→0) log N(δ) / log(1/δ)。实际操作中有四点要注意第一目标需要二值化取值只有0和1第二δ的序列通常取2的幂次例如2、4、8……直到接近图像尺寸第三统计N(δ)时只要盒子内有任意一个非零像素就计数第四在双对数坐标下用最小二乘拟合斜率而不是用最后两点直接连线因为直接连线对边界像素和尺度选择的噪声都非常敏感。实际计算时我会剔除尺寸过大和过小的尺度段过小的尺度会被像素分辨率限制过大的尺度则没有足够盒子样本拟合区间通常取δ2到图像短边的一半。2.2 多重分形谱的数学框架配分函数与Legendre变换多重分形分析的第一步是把目标转换为概率测度。假设p_i(δ)表示第i个盒子内的质量占总质量的比重定义配分函数Z(q,δ)Σ p_i(δ)^q。当δ趋于0时Z(q,δ)以幂律形式δ^τ(q)变化这里的τ(q)称为质量指数通过对log Z(q,δ)与log δ拟合得到。在得到τ(q)后通过Legendre变换得到奇异谱α(q)dτ(q)/dqf(q)q·α(q)-τ(q)。这里所有公式都只依赖于τ(q)的斜率。此外还能得到广义维数D(q)当q0时D(0)-τ(0)就是容量维数也就是盒子维数q1时D(1)对应信息维数q2时D(2)对应关联维数。需要特别留意D(1)不能直接用τ(1)/(1-1)算因为q1处τ(1)0且分母也为0要用α(1)来替代。2.3 q值到底有什么用权重因子的物理含义很多人在代码里看到q从-10取到10不知道为什么要这么干。q本质上是一个权重因子当q0时p_i^q会放大概率较大的盒子反映测度集中、高密度区域的行为当q0时p_i^q会放大概率接近0的盒子反映测度稀疏、奇异性较强的局部区域。因此要完整刻画多重分形q必须在负数和正数两侧都取到。这里的符号约定在不同文献里略有差别我的习惯是先按Z≈δ^τ的规则拟合得到αf(q)再算f(q)q·α-τ最后用已知测度做一次基准验证确保曲线形状和理论一致否则很可能是符号或负号写反了。3. MATLAB工程实现可直接复用的核心代码3.1 盒子维数计算函数实现先给出直接能抄的盒计数法函数。输入是二值图像输出是盒子维数D和用于调试的坐标数据。function [D, logInvScale, logCount] boxcount(bwImg) bwImg logical(bwImg); [Nx, Ny] size(bwImg); maxLevel floor(log2(max(Nx, Ny))); scales 2.^(1:maxLevel); nS numel(scales); logInvScale zeros(nS, 1); logCount zeros(nS, 1); for k 1:nS s scales(k); nBoxX ceil(Nx / s); nBoxY ceil(Ny / s); padImg false(nBoxX * s, nBoxY * s); padImg(1:Nx, 1:Ny) bwImg; gridCell reshape(padImg, s, nBoxX, s, nBoxY); hit squeeze(any(any(gridCell, 1), 3)); logInvScale(k) log(1 / s); logCount(k) log(sum(hit(:))); end valid isfinite(logCount) (logCount -inf); p polyfit(logInvScale(valid), logCount(valid), 1); D p(1); end这段代码按网格切分图像并统计每个盒子内是否有目标点polyfit做一阶拟合得到斜率。MATLAB的log默认是自然对数不影响斜率结果。图像尺寸不是盒子边长的整数倍时代码用零填充补齐空盒不会干扰计数。实际使用时可以在尺度序列两侧各截掉一层规避边界效应。3.2 多重分形谱算法的完整代码多重分形谱的代码化并不复杂但每一步都容易出细节错误。把核心流程展开function [alpha, falpha, Dq, qRange] multifractal_spectrum(bwImg, qRange, nScales) img double(bwImg); [Nx, Ny] size(img); tot sum(img(:)); if tot 0 error(输入图像全黑无法计算); end p img / tot; scaleList unique(round(logspace(log10(2), log10(min(Nx, Ny) / 4), nScales))); nSc numel(scaleList); nQ numel(qRange); logZ zeros(nSc, nQ); for si 1:nSc s scaleList(si); nBoxX ceil(Nx / s); nBoxY ceil(Ny / s); padImg zeros(nBoxX * s, nBoxY * s); padImg(1:Nx, 1:Ny) p; gridCell reshape(padImg, s, nBoxX, s, nBoxY); boxMass squeeze(sum(sum(gridCell, 1), 3)); boxMass boxMass(:); boxMass boxMass(boxMass 0); for qi 1:nQ logZ(si, qi) log(sum(boxMass.^qRange(qi))); end end logEps log(1 ./ scaleList); tau zeros(nQ, 1); alpha zeros(nQ, 1); falpha zeros(nQ, 1); Dq zeros(nQ, 1); for qi 1:nQ pFit polyfit(logEps, logZ(:, qi), 1); tau(qi) pFit(1); end alpha gradient(tau) ./ gradient(qRange(:)); falpha qRange(:) .* alpha - tau; for qi 1:nQ if abs(qRange(qi) - 1) 1e-8 Dq(qi) alpha(qi); else Dq(qi) tau(qi) / (qRange(qi) - 1); end end end这里用logspace生成尺度序列尺度上限定为图像短边的四分之一避免网格过少。对每个尺度统计每个盒子里的概率质量空盒直接排除这是避免log(0)的关键。拟合时回归变量是log(1/scale)对应前文约定的符号规则。输出的alpha、falpha就是多重分形谱的横纵坐标Dq是广义维数谱。3.3 参数选型尺度范围、q范围与拟合区间参数选型的学问比代码本身更大。尺度范围如果上限太大盒子数量太少统计噪声很严重上限太小又捕捉不到大尺度上的标度行为。q范围建议至少取[-10, 10]如果目标概率测度分布特别不均匀再扩展到[-20, 20]否则谱的两端会明显缺失。拟合时不要使用全部尺度点可以观察双对数散点图去掉两端明显偏离直线的点。在调试时我会先输出logEps和logZ用plot看一下线性回归的残差分布而不是盲信polyfit的拟合优度。4. 验证与结果解读代码跑通只是第一步4.1 用Sierpinski三角形和Cantor集验证算法代码写完后第一件事是验证。我最常用的是Cantor集和Sierpinski三角形这两个对象有解析理论值。Sierpinski三角形的分形维数是log3/log2≈1.5850我用512×512的二值图跑盒计数法多次实验都在1.55~1.60之间与理论值误差在2%以内。这个误差主要来源于有限分辨率和边界像素属于正常情况。如果误差超过5%就要检查二值化是否正确、尺度序列是否取到了边界效应明显的层。用Cantor集验证多重分形谱时由于它是均匀分形f(α)曲线会非常窄近似集中于理论维数log2/log3≈0.6309附近如果计算出的曲线明显变宽就要怀疑配分函数或拟合符号出了问题。4.2 f(α)曲线和Dq谱到底该怎么读通常f(α)是一条上凸的钟形曲线峰值对应的纵坐标就是容量维数D0即盒子维数。曲线的左右端点分别对应q趋于正无穷和负无穷的行为所以当q范围不够大时曲线两端会显得“缺角”。曲线越宽说明对象的局部奇异性差异越大多重分形特征越强曲线越窄说明对象越接近均匀分形。Dq谱随着q增大单调不增特殊点D0D1D2是正常现象。如果算出的Dq谱出现上升趋势或剧烈震荡基本可以断定拟合尺度范围选取不当或者图像噪声干扰过大。4.3 实际数据案例分析从曲线形态看结构差异以岩石薄片孔隙图像为例我曾算过一组数据D0约1.72D1约1.64D2约1.55谱宽也明显大于0说明孔隙分布具有显著多重分形特征局部高密度区域和低密度区域并存。而用同一套代码处理一幅规则网格图像时D0、D1、D2几乎相等f(α)曲线缩成一个很窄的峰判断为均匀结构。这个对比非常直观实际分析报告里我一般同时给出D0、D1、D2和谱宽四个量而不是只给一条f(α)曲线这样更便于不同样本之间的横向比较。5. 实际场景、避坑指南与实操心得5.1 这套算法能用在哪些地方多重分形谱和盒子维数在工业界用得比想象中广。图像纹理分析方面医学影像中的组织病变区域、岩石铸体薄片中的孔隙结构都能通过f(α)曲线的宽度和峰值来判断复杂度变化信号处理方面机械振动信号的多重分形谱可以用于故障特征提取磨损、裂纹等异常状态往往表现为Dq谱的明显变化此外还有金融时间序列分析把收益率序列的波动转化为概率测度后做多重分形分析能刻画市场在不同时间尺度上的波动聚集性。只要能把目标数据转换为非负的测度比如灰度、能量、质量、频率这套算法就可以迁移过去。5.2 常见问题排查一份速查表在实际跑代码时我遇到过反馈较多的四类问题整理成速查表问题可能原因处理方法盒子维数结果偏离理论值二值化不准确或边界效应检查图像预处理去掉过大过小尺度再拟合log中出现NaN或Inf盒子概率为0时取log统计前先排除零盒子或对像素值加极小epsilonf(α)谱两端缺失严重q范围太小或尺度区间不足扩大q范围增加尺度序列密度Dq谱震荡不单调拟合尺度区间包含非线性段手动调整拟合区间用残差图指导选取这里需要特别强调加epsilon的做法要谨慎因为epsilon太大会改变概率测度的分布特性从而污染谱的真实形状我一般优先选择排除空盒的方式而不是加噪声。5.3 几条实操心得最后分享一些我的个人习惯。第一不要迷信单次计算结果分形维数对尺度范围非常敏感报告结果时必须附带拟合区间和散点图很多审稿人会追问这一点。第二多重分形谱的拟合建议用最小二乘前先画图人眼判断比任何统计指标都可靠。第三数据量小时不要强行算多重分形至少要保证最小尺度上还有足够的盒子数否则会出现严重的统计波动一般图像尺寸不要小于256×256信号序列则建议长度大于1024。第四如果只是需要一个大致的复杂度量化指标盒子维数就够了只有当你需要了解局部非均匀性时再上多重分形谱避免杀鸡用牛刀。本文还有配套的精品资源点击获取