简介：这份资源面向分形几何与复杂系统分析的学习者与研究者，提供Matlab环境下的多重分形谱算法与盒子维数计算实现，适合具备一定数值计算基础、希望动手复现分形维数估计的读者。压缩包共2个文件，均为.m脚本，整体约3KB，分别对应多重分形谱计算与盒子分形维数计算两条主线，可直接在Matlab中运行与修改。资源围绕分箱细化、盒计数统计、H(q)与f(α)特征谱构造，以及log-log曲线拟合求盒子维数等关键环节展开，帮助读者理解从数据准备到维数输出的完整流程。目前已有1343人学习下载，可作为图像处理、信号分析、金融建模等场景中分形特征提取的入门参考与代码模板，便于快速搭建实验并对照验证结果。

1. 从一组“长得不像规则图形”的数据说起：这套多重分形谱+盒子维数代码到底能干什么

如果你手头有一批表面形貌、纹理图像、时间序列或者空间分布数据，领导让你“量化一下它的复杂程度”，你大概率会先想到分形维数。但单一的分形维数只能告诉你“有多复杂”，没法告诉你“复杂得均不均匀”。这时候多重分形谱就派上用场了——它能把一个分形对象拆成不同奇异强度的子集，给出一个谱函数 f(α)，横轴 α 是奇异指数，纵轴 f(α) 是分形维数。而盒子维数则是计算多重分形谱的基础工具，负责在不同尺度下统计覆盖盒子数。

这套 MATLAB 代码包解决的就是这件事：输入一张二值图或灰度图，输出盒子维数曲线和多重分形谱。适合做材料表面分析、生物医学图像纹理、地质孔隙结构、信号复杂度评估的从业者。代码不依赖额外工具箱，纯 MATLAB 实现，拿到就能跑。下面我从原理到实操拆一遍，顺带把几个容易翻车的地方讲清楚。

2. 盒子维数与多重分形谱：原理、参数和代码结构

2.1 盒子维数的计算逻辑与尺度选择

盒子维数的定义很直白：用边长为 ε 的盒子去覆盖目标集合，统计非空盒子数 N(ε)，然后看 N(ε) 随 ε 减小怎么增长。如果满足 N(ε) ∝ ε^(-D)，那 D 就是盒子维数。实际计算时，对 log N(ε) 和 log(1/ε) 做线性拟合，斜率就是 D。

但这里有个关键参数：ε 的取值序列。常见做法是取 2 的幂次，比如 ε = 2, 4, 8, 16, 32, 64，这样盒子网格能整齐划分图像。如果图像尺寸不是 2 的幂，要么补零，要么取最大公约数序列。我一般会先检查图像尺寸，不是 2 的幂就裁剪或填充到最近的 2 的幂，避免网格偏移导致统计误差。

代码里对应的核心函数是boxdim.m，输入是二值矩阵，输出是拟合斜率。下面这段是尺度序列生成和盒子计数的核心逻辑：

% 输入：BW 为二值图像矩阵，要求尺寸为 2 的幂 % 输出：D 为盒子维数，logN 和 logInvEps 用于绘图验证 function [D, logN, logInvEps] = boxdim(BW) [h, w] = size(BW); maxPow = log2(min(h, w)); epsList = 2.^(1:floor(maxPow)); % 尺度序列：2,4,8,... N = zeros(size(epsList)); for k = 1:length(epsList) eps = epsList(k); count = 0; for i = 1:eps:h for j = 1:eps:w block = BW(i:min(i+eps-1,h), j:min(j+eps-1,w)); if any(block(:)) count = count + 1; end end end N(k) = count; end logN = log(N); logInvEps = log(1./epsList); p = polyfit(logInvEps, logN, 1); D = p(1); end

逻辑说明：外层循环遍历尺度，内层双重循环按 eps 步长切块，any(block(:))判断块内是否有前景像素。参数方面，epsList从 2 开始而不是 1，因为 eps=1 时每个像素一个盒子，N 等于前景像素总数，对拟合贡献不大且容易受噪声影响。拟合用一次多项式，斜率即 D。如果拟合的 R² 低于 0.98，说明尺度范围选得不好，需要调整 epsList 的上下限。

2.2 多重分形谱的配分函数与 Legendre 变换

多重分形谱的核心是配分函数 χ(q, ε) = Σ μ_i^q，其中 μ_i 是第 i 个盒子的测度（比如像素灰度占比），q 是矩阶数。对每个 q，计算 ln χ(q, ε) 对 ln ε 的斜率，得到质量指数 τ(q)。然后通过 Legendre 变换：α = dτ/dq，f(α) = qα - τ(q)，得到多重分形谱。

q 的取值范围很关键。q 太大或太小都会导致数值溢出或拟合不稳定。常见做法是 q 从 -10 到 10，步长 0.5，共 41 个点。如果数据动态范围大，可以扩展到 -20 到 20，但要注意 μ_i^q 可能下溢或上溢，需要做对数域处理。

代码里mfspectrum.m负责这部分，核心片段：

% 输入：BW 二值图或灰度图（归一化到 0-1） % 输出：alpha 和 falpha 为多重分形谱坐标 function [alpha, falpha, qList] = mfspectrum(img) img = double(img); img = img / sum(img(:)); % 归一化为概率测度 [h, w] = size(img); maxPow = log2(min(h, w)); epsList = 2.^(2:floor(maxPow)); % 尺度从 4 开始，避免过小尺度噪声 qList = -10:0.5:10; tau = zeros(size(qList)); for qi = 1:length(qList) q = qList(qi); lnChi = zeros(size(epsList)); for k = 1:length(epsList) eps = epsList(k); chi = 0; for i = 1:eps:h for j = 1:eps:w block = img(i:min(i+eps-1,h), j:min(j+eps-1,w)); mu = sum(block(:)); if mu > 0 chi = chi + mu^q; end end end lnChi(k) = log(chi); end p = polyfit(log(1./epsList), lnChi, 1); tau(qi) = p(1); end % Legendre 变换 alpha = diff(tau) ./ diff(qList); falpha = qList(1:end-1) .* alpha - tau(1:end-1); end

逻辑说明：先归一化图像使总测度为 1，保证 μ_i 是概率。对每个 q，计算不同尺度下的配分函数，取对数后对 log(1/ε) 拟合斜率得 τ(q)。最后用差分近似导数得到 α 和 f(α)。参数注意：epsList从 4 开始，因为尺度太小盒子内测度接近 0，log 会出问题；q 范围根据数据调整，如果谱曲线两端翘起明显，说明 q 范围不够。

2.3 代码包的文件结构与调用关系

这套代码包通常包含这几个文件：boxdim.m算盒子维数，mfspectrum.m算多重分形谱，main_demo.m是示例脚本，sample_data.mat或示例图片。调用顺序是先用boxdim验证数据的分形特性，再用mfspectrum看谱形。如果盒子维数拟合 R² 很低，多重分形谱的结果也不可信。

我一般会在main_demo.m里加一段可视化，把 logN-logInvEps 散点和拟合线画出来，肉眼确认线性关系。谱曲线则画 α-f(α) 图，正常的多重分形谱应该是单峰钟形，如果出现多峰或断裂，多半是 q 范围或尺度范围有问题。

3. 实操：从图像预处理到谱曲线输出的完整流程

3.1 图像预处理与二值化阈值选择

拿到一张原始图像，不能直接扔进mfspectrum。先要转灰度、去噪、二值化。二值化阈值的选择直接影响盒子计数和测度分布。常见做法是用 Otsu 方法自动阈值，或者根据直方图手动选。如果图像背景不均匀，先做背景校正或局部阈值。

% 读取图像并预处理 raw = imread('sample.png'); if size(raw, 3) == 3 gray = rgb2gray(raw); else gray = raw; end gray = im2double(gray); % 中值滤波去椒盐噪声 gray = medfilt2(gray, [3 3]); % Otsu 二值化 level = graythresh(gray); BW = imbinarize(gray, level); % 裁剪到 2 的幂尺寸 [h, w] = size(BW); target = 2^floor(log2(min(h, w))); BW = BW(1:target, 1:target);

逻辑说明：rgb2gray处理彩色图，medfilt2去噪避免孤立像素影响盒子计数，graythresh返回 Otsu 阈值，imbinarize生成二值图。最后裁剪到 2 的幂，保证网格整齐。参数注意：中值滤波窗口用 3×3 就够，太大糊掉细节；如果图像本身是二值的，跳过二值化直接裁剪。

3.2 盒子维数计算与拟合质量检查

预处理完先跑boxdim，看拟合质量。如果 R² 低于 0.98，别急着往下走，先检查尺度序列和图像质量。

[D, logN, logInvEps] = boxdim(BW); % 拟合质量检查 p = polyfit(logInvEps, logN, 1); yfit = polyval(p, logInvEps); SSres = sum((logN - yfit).^2); SStot = sum((logN - mean(logN)).^2); R2 = 1 - SSres / SStot; fprintf('盒子维数 D = %.4f, R2 = %.4f\n', D, R2); % 绘图验证 figure; plot(logInvEps, logN, 'o', logInvEps, yfit, '-'); xlabel('log(1/\epsilon)'); ylabel('log N(\epsilon)'); title(sprintf('Box Dimension Fit, D=%.4f, R^2=%.4f', D, R2));

逻辑说明：polyfit做一次拟合，手动算 R² 判断线性度。绘图是为了肉眼确认散点是否落在直线附近。如果 R² 低，常见原因是图像噪声大、尺度范围太窄或太宽、或者数据本身不是分形。我一般会尝试去掉最小和最大尺度各一个点再拟合，往往能改善。

3.3 多重分形谱计算与谱曲线解读

盒子维数过关后，跑mfspectrum。注意输入应该是归一化后的灰度图或二值图，不是原始图像。

[alpha, falpha, qList] = mfspectrum(BW); figure; plot(alpha, falpha, 'b-', 'LineWidth', 1.5); xlabel('\alpha'); ylabel('f(\alpha)'); title('Multifractal Spectrum'); grid on; % 谱宽和对称性指标 spectrumWidth = max(alpha) - min(alpha); [~, idx] = max(falpha); asymmetry = (alpha(idx) - min(alpha)) / (max(alpha) - alpha(idx)); fprintf('谱宽 = %.4f, 不对称度 = %.4f\n', spectrumWidth, asymmetry);

逻辑说明：alpha和falpha是谱坐标，qList是矩阶序列。谱宽反映复杂度的动态范围，越宽说明奇异性分布越不均匀。不对称度反映谱形偏向哪一侧，左偏表示低奇异指数主导，右偏表示高奇异指数主导。参数注意：如果谱曲线出现明显截断或尖峰，检查 q 范围是否覆盖了 τ(q) 的线性区间。

4. 避坑与排查：盒子计数和谱拟合里最容易翻车的五个点

4.1 现象：盒子维数拟合 R² 只有 0.8 左右，散点明显弯曲

原因：尺度序列选得太宽，小尺度受噪声主导，大尺度盒子数太少统计不稳。或者图像本身不是分形，强行拟合没有意义。

解决：把epsList限制在中间段，比如从 4 到 64，去掉最小和最大尺度。同时检查图像预处理是否过度滤波导致细节丢失。如果 R² 还是低，换一组数据验证代码本身没问题。

4.2 现象：多重分形谱出现负的 f(α) 值

原因：Legendre 变换用差分近似导数时，q 步长太大或 τ(q) 拟合误差被放大。另外配分函数计算时 μ_i 有零值，log 处理不当。

解决：把 q 步长从 0.5 改成 0.2，增加 q 点数。配分函数里对 μ_i=0 的盒子跳过，不要加 log(0)。如果还有负值，检查 τ(q) 是否单调递减，正常应该单调。

4.3 现象：谱曲线只有半个峰，右边或左边被截断

原因：q 范围不够。q 正得不够大，高奇异指数部分没覆盖；q 负得不够小，低奇异指数部分没覆盖。

解决：把 q 范围扩展到 -20 到 20，步长 0.5。但要注意数值溢出，μ_i^q 在 q 很大时可能下溢到 0，需要在对数域计算：用q * log(mu)代替log(mu^q)。

4.4 现象：同一张图跑两次结果不一样

原因：图像读取时格式不同导致灰度范围不同，或者二值化阈值随机。另外如果代码里有随机数（比如随机采样），也会导致结果波动。

解决：固定图像读取后的归一化方式，二值化用确定阈值。检查代码里有没有rand或randn，有的话设随机种子rng(0)。我一般会在脚本开头加rng(0)，保证可复现。

4.5 现象：大图像跑得特别慢，内存爆掉

原因：双重循环按像素遍历，图像尺寸 2048×2048 时循环次数太多。另外配分函数对每个 q 都重新切块，重复计算。

解决：把盒子切块索引预先算好，或者用blockproc替代手动循环。更彻底的做法是向量化：用accumarray按盒子编号累加测度。如果内存不够，先把图像降采样到 1024×1024 再跑，分形维数对降采样不敏感。

5. 进阶技巧：用滑动窗口做局部多重分形谱，验证代码鲁棒性

全局多重分形谱给的是整张图的平均统计，但很多场景下我们想知道不同区域的复杂度差异。比如材料表面有的区域粗糙有的区域平滑，全局谱会把这些差异平均掉。这时候可以用滑动窗口，把图像切成小块，每块单独算谱，然后看谱宽和不对称度的空间分布。

具体做法：设定窗口大小比如 256×256，步长 128，对每个窗口跑mfspectrum，记录谱宽和 α 的均值。最后画一张热力图，横纵坐标是窗口位置，颜色是谱宽。这样能直观看到哪些区域复杂度高。

% 滑动窗口多重分形谱 winSize = 256; step = 128; [h, w] = size(BW); rows = 1:step:(h - winSize + 1); cols = 1:step:(w - winSize + 1); spectrumMap = zeros(length(rows), length(cols)); for ri = 1:length(rows) for ci = 1:length(cols) block = BW(rows(ri):rows(ri)+winSize-1, cols(ci):cols(ci)+winSize-1); [alpha, ~, ~] = mfspectrum(block); spectrumMap(ri, ci) = max(alpha) - min(alpha); end end figure; imagesc(spectrumMap); colorbar; title('Local Multifractal Spectrum Width'); xlabel('Window Column'); ylabel('Window Row');

逻辑说明：外层双重循环遍历窗口位置，每个窗口调用mfspectrum，只取谱宽存入矩阵。imagesc画热力图，颜色越亮谱宽越大。参数注意：窗口大小要是 2 的幂，否则mfspectrum里的尺度序列会出问题。步长一般取窗口一半，保证重叠。如果窗口太小比如 64×64，尺度序列只有 2 和 4 两个点，拟合不可靠，所以最小窗口建议 128。

验证代码鲁棒性的方法：拿一张已知分形维数的合成图，比如 Sierpinski 三角形或 Cantor 集，跑一遍看结果是否接近理论值。Sierpinski 三角形的盒子维数是 log3/log2 ≈ 1.585，多重分形谱应该是一个点（因为它是单分形）。如果代码跑出来谱宽接近 0，说明实现正确。如果谱宽很大，说明配分函数或 Legendre 变换有问题。

我自己的习惯是：每次拿到新的分形计算代码，先用合成图验证，再用真实数据跑。合成图验证通过后，真实数据的结果才有可信度。从那以后我每次跑多重分形谱之前，都强制走一遍合成图验证流程，省得后面返工。希望帮到你。

本文还有配套的精品资源，点击获取