简介：多重分形谱算法与盒子维数计算是分形几何中分析复杂系统自相似结构的重要工具，这套基于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)/dq，f(q)=q·α(q)-τ(q)。这里所有公式都只依赖于τ(q)的斜率。此外还能得到广义维数D(q)，当q=0时D(0)=-τ(0)就是容量维数，也就是盒子维数；q=1时D(1)对应信息维数；q=2时D(2)对应关联维数。需要特别留意，D(1)不能直接用τ(1)/(1-1)算，因为q=1处τ(1)=0且分母也为0，要用α(1)来替代。
2.3 q值到底有什么用：权重因子的物理含义
很多人在代码里看到q从-10取到10，不知道为什么要这么干。q本质上是一个权重因子：当q>0时，p_i^q会放大概率较大的盒子，反映测度集中、高密度区域的行为；当q<0时，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增大单调不增，特殊点D0>D1>D2是正常现象。如果算出的Dq谱出现上升趋势或剧烈震荡，基本可以断定拟合尺度范围选取不当，或者图像噪声干扰过大。
4.3 实际数据案例分析：从曲线形态看结构差异
以岩石薄片孔隙图像为例，我曾算过一组数据：D0约1.72，D1约1.64，D2约1.55，谱宽也明显大于0，说明孔隙分布具有显著多重分形特征，局部高密度区域和低密度区域并存。而用同一套代码处理一幅规则网格图像时，D0、D1、D2几乎相等，f(α)曲线缩成一个很窄的峰，判断为均匀结构。这个对比非常直观，实际分析报告里我一般同时给出D0、D1、D2和谱宽四个量，而不是只给一条f(α)曲线，这样更便于不同样本之间的横向比较。
5. 实际场景、避坑指南与实操心得
5.1 这套算法能用在哪些地方
多重分形谱和盒子维数在工业界用得比想象中广。图像纹理分析方面，医学影像中的组织病变区域、岩石铸体薄片中的孔隙结构，都能通过f(α)曲线的宽度和峰值来判断复杂度变化；信号处理方面，机械振动信号的多重分形谱可以用于故障特征提取，磨损、裂纹等异常状态往往表现为Dq谱的明显变化；此外还有金融时间序列分析，把收益率序列的波动转化为概率测度后做多重分形分析，能刻画市场在不同时间尺度上的波动聚集性。只要能把目标数据转换为非负的测度，比如灰度、能量、质量、频率，这套算法就可以迁移过去。
5.2 常见问题排查：一份速查表
在实际跑代码时，我遇到过反馈较多的四类问题，整理成速查表：
| 问题 | 可能原因 | 处理方法 |
|---|---|---|
| 盒子维数结果偏离理论值 | 二值化不准确或边界效应 | 检查图像预处理，去掉过大过小尺度再拟合 |
| log中出现NaN或Inf | 盒子概率为0时取log | 统计前先排除零盒子，或对像素值加极小epsilon |
| f(α)谱两端缺失严重 | q范围太小或尺度区间不足 | 扩大q范围，增加尺度序列密度 |
| Dq谱震荡不单调 | 拟合尺度区间包含非线性段 | 手动调整拟合区间，用残差图指导选取 |
这里需要特别强调，加epsilon的做法要谨慎，因为epsilon太大会改变概率测度的分布特性，从而污染谱的真实形状；我一般优先选择排除空盒的方式，而不是加噪声。
5.3 几条实操心得
最后分享一些我的个人习惯。第一，不要迷信单次计算结果，分形维数对尺度范围非常敏感，报告结果时必须附带拟合区间和散点图，很多审稿人会追问这一点。第二，多重分形谱的拟合建议用最小二乘前先画图，人眼判断比任何统计指标都可靠。第三，数据量小时不要强行算多重分形，至少要保证最小尺度上还有足够的盒子数，否则会出现严重的统计波动；一般图像尺寸不要小于256×256，信号序列则建议长度大于1024。第四，如果只是需要一个大致的复杂度量化指标，盒子维数就够了；只有当你需要了解局部非均匀性时，再上多重分形谱，避免杀鸡用牛刀。
本文还有配套的精品资源，点击获取
