简介：分形维数计算是复杂系统与非规则形状研究中的重要工具，在图像处理、信号分析、材料形貌表征等领域应用广泛。这份面向科研工作者与学习者的代码包，集中实现了差分盒维数、功率谱分析及结构函数三类主流算法，可帮助研究者快速估算不规则对象的复杂程度。资源共包含五个m脚本文件，压缩包整体仅约3KB，各脚本按算法模块独立组织，结构清晰易读，便于对照学习与二次开发，目前已有近三千人次学习下载。运行代码后，读者可掌握差分盒维数通过盒子统计与对数拟合估算分形维数的思路，借助功率谱分析探索信号的频率成分与统计特性，利用结构函数考察空间数据在不同尺度下的变化规律，进而形成一套完整的分形量化分析流程，适合正在学习分形理论并需要快速开展数值实验的研究者参考。
1. 分形维数计算：为什么不同算法会算出不同结果
做材料断面、表面粗糙度或图像纹理分析时，一旦涉及“这个表面到底有多不规则”，绕不开的就是分形维数。分形维数计算的matlab实现里，三个词出现频率最高：差分盒维数、功率谱、结构函数。它们都能给“不规则程度”一个数值，但同一个样本，用三种方法算出来的结果可能差 0.2 甚至 0.4，这不是代码 bug，而是各自对“几何尺度”的定义不同。我的经验是：理解这一点比会跑代码更重要。
这篇文章面向两类读者：一类是刚接触分形理论、想用 Matlab 算出第一张图像的维数；另一类是做粗糙度或纹理分析、已经跑出结果但不知道怎么调参数、怎么解释数值差异。全文从差分盒维数讲起，再做功率谱和结构函数，最后给一套验证和避坑清单。看完你能独立跑通三个方向，并且能判断哪个结果可信、哪个只是巧合。
2. 差分盒维数：Matlab 实现的核心算法
2.1 为什么差分盒维数成为默认选择
在图像纹理和表面形貌分析里，差分盒维数（Differential Box Counting，DBC）几乎是默认首选。原因很实际：它对灰度图像直接操作，不需要额外提取轮廓或二值化，把一个二维灰度曲面当成三维几何体来数盒子。实现逻辑和经典盒计数一致，但差分之处在于：它用一个动态的盒子高度去覆盖灰度方向的变化，而不是固定高度。对于表面粗糙度分析，灰度值的起伏恰好对应表面高度变化，这种“动态高度”的做法非常贴合数据形态。
选 DBC 还有个重要理由——抗噪能力相对好。它统计的是不同尺度下盒子数量的对数斜率，局部噪声被淹没在统计过程里。相比逐像素梯度类方法，DBC 对轻微的像素抖动不敏感，这对 SEM 图像、金相图像、磨损表面这类真实样本尤其重要。在我接触过的方案里，十有八九是先用 DBC 算出一个全局维数，再用其他方法做交叉验证。
2.2 差分盒维数的 Matlab 核心代码
我先给出一份可以直接落地的函数实现。这个实现遵循经典 DBC 流程：灰度归一化、确定盒子尺寸序列、统计各尺度下盒子数量、对数域线性拟合。
function D = dbc_fractal_dim(I, box_sizes) % 差分盒维数计算 % I : 灰度图像, double 类型, 范围 [0, 1] % box_sizes : 盒子边长序列, 例如 [2 4 8 16 32 64] % D : 估计的分形维数 I = im2double(I); [rows, cols] = size(I); L = min(rows, cols); % 有效边长, 保证网格为正方形 max_gray = max(I(:)); min_gray = min(I(:)); G = 256; % 灰度级数, 归一化到 [0,1] 后映射到 [1,G] n_g = round((G / (max_gray - min_gray)) * (I - min_gray)); n_g = max(n_g, 1); % 避免除零和灰度为0的情况 Nr = zeros(size(box_sizes)); % 盒子计数结果 for k = 1:length(box_sizes) s = box_sizes(k); if s > L/2 continue; % 盒子边长超过图像一半, 统计失去意义 end grid_r = floor(rows / s); grid_c = floor(cols / s); h = floor(max_gray * s / G) + 1; % 盒子高度随灰度动态调整 count = 0; for i = 1:grid_r for j = 1:grid_c % 当前网格灰度最小值和最大值 block = n_g((i-1)*s+1 : i*s, (j-1)*s+1 : j*s); min_val = min(block(:)); max_val = max(block(:)); count = count + (floor(max_val / h) - floor(min_val / h) + 1); end end Nr(k) = count; end % 只用有效尺寸做拟合 valid = (Nr > 0); x = log(1 ./ box_sizes(valid)); y = log(Nr(valid)); p = polyfit(x, y, 1); D = p(1); end这段代码里最值得说明的是h = floor(max_gray * s / G) + 1这一行。差分盒维数的核心就在于：当盒子边长s增大时，灰度方向的盒子高度也按比例放大，这样才能保持盒子在三维空间中的“立方体”属性。如果把高度固定不变，那么小盒子会数出过多的盒子数，大盒子又数得过少，拟合斜率会偏大，结果虚高。+1是为了防止灰度为 0 时出现零高度。
盒子尺寸序列box_sizes是另一个关键参数。常见做法是取 2 的整数次幂，例如[2 4 8 16 32 64]。但要注意，s必须大于等于 2，因为边长为 1 时每个盒子只覆盖一个像素，盒子数等于像素数，维数恒等于 2，没有任何区分度。另外，s超过图像边长的一半后，网格数太少，统计意义丧失，代码里已经用continue跳过。
2.3 从二维到一维：让 DBC 适配轮廓或剖面线
图像整体维数之外，很多场景需要算某条剖面线的维数。比如表面粗糙度分析里，一条轮廓线的分形维数比整幅图的维数更能反映加工工艺的细节。DBC 同样可以用于一维信号，只需把二维网格退化成一维区间。
function D = dbc_1d(signal, box_sizes) % 一维信号差分盒维数 % signal: 行向量, 代表表面轮廓高度或灰度剖面 % box_sizes: 区间长度序列 signal = signal(:)'; N = length(signal); max_val = max(signal); min_val = min(signal); G = 256; s_norm = round((G / (max_val - min_val)) * (signal - min_val)); Nr = zeros(size(box_sizes)); for k = 1:length(box_sizes) s = box_sizes(k); segs = floor(N / s); h = floor(max_val * s / G) + 1; count = 0; for i = 1:segs seg = s_norm((i-1)*s+1 : i*s); count = count + (floor(max(seg) / h) - floor(min(seg) / h) + 1); end Nr(k) = count; end valid = Nr > 0; x = log(1 ./ box_sizes(valid)); y = log(Nr(valid)); D = polyfit(x, y, 1); D = D(1); end一维版本的原理和二维完全一致，只是把“块”换成“段”。这里有个容易翻车的地方：一维信号的维数理论范围是 1 到 2，但 DBC 结果经常略高于 2，尤其是信号长度只有几百个点的时候。原因是边界上的离群点把小尺度下的盒子数抬高，拟合斜率随之变大。处理办法是在拟合时去掉最小的 1 到 2 个尺寸，或者对信号做轻度平滑。
3. 功率谱与结构函数法：从频域和时域两个方向逼近
3.1 功率谱法：从频域斜率反推维数
差分盒维数在空间域数盒子，功率谱法走完全不同的路径。它把信号或图像变换到频域，利用分形信号功率谱的幂律特性：功率谱密度P(f)与频率f在对数坐标下呈线性关系，斜率 β 与分形维数 D 之间存在解析关系。对于一维表面轮廓，D = (5 - β) / 2；对于二维表面，D = (8 - β) / 2（其中 β 是径向平均功率谱的斜率）。这个方法的优势在于：它能区分不同方向上的分形特征，并且对周期分量非常敏感——如果信号里混有周期性纹理，功率谱会出现尖峰，一眼就能看出来。
Matlab 实现功率谱法的核心是fft2配合径向平均。下面给出一维信号和二维图像两个版本，实际使用中一维版本更常用于表面轮廓仪数据。
function D = psd_fractal_dim_1d(signal, fs) % 功率谱法计算一维信号分形维数 % signal: 表面轮廓高度, 行向量 % fs : 采样率, 用于频率轴, 不参与维数计算但便于截断 signal = signal(:)'; N = length(signal); % 去趋势: 避免直流分量和线性趋势主导低频 x = (0:N-1)'; p = polyfit(x, signal', 1); signal_detrend = signal' - polyval(p, x); % Welch 法估计功率谱, 分段平均能降低方差 [psd, freq] = pwelch(signal_detrend, hann(N/4), N/8, N, fs); % 截取有效频段: 去掉直流和极高频率 f_min = fs / N * 4; % 至少覆盖 4 个周期 f_max = fs / 2 * 0.8; % 避开奈奎斯特附近的衰减 idx = freq >= f_min & freq <= f_max; x = log10(freq(idx)); y = log10(psd(idx)); p = polyfit(x, y, 1); beta = -p(1); % 功率谱斜率, 取负号 D = (5 - beta) / 2; end功率谱法的参数陷阱全在频段截取上。如果f_min取得太小，低频段的功率谱受趋势项影响，斜率会明显偏小；如果f_max取得太大，高频段进入噪声平台，斜率趋近于 0，算出的维数直接跌穿理论下限。我一般会画一次loglog(freq, psd)看曲线：选择斜率稳定的中间段做拟合，而不是盲目套固定比例。上面代码里f_min = fs/N*4和f_max = fs/2*0.8只是保守起点，真实数据十有八九需要手动调整。
二维图像版本的功率谱法，处理流程是：对图像加窗（比如 Hann 窗）消除边缘不连续、做fft2、计算功率谱、按径向距离做环形平均，最后在对数坐标拟合。要注意的是，二维表面如果各向异性明显，径向平均会掩盖方向特征，这时候应该把角度分成几个扇区分别拟合，得到的多个维数可以用于评估表面各向异性程度。这也是功率谱法对比 DBC 的一个核心优势：DBC 只能给一个全局数值，而功率谱法能给出方向维数。
3.2 结构函数法：从高度差统计中提取 Hurst 指数
结构函数法是我个人在表面粗糙度分析里用得最多的方法，原因是它直接在空间域操作，对非平稳信号更稳健。结构函数定义为：S(τ) = E[(z(x+τ) - z(x))^2]，它描述的是相隔距离 τ 的两点高度差的均方值。对于分形表面，S(τ)与 τ 在对数坐标下线性相关，斜率为2H，其中 H 是 Hurst 指数，分形维数D = 2 - H（一维轮廓）或D = 3 - H（二维表面）。
function [D, H] = structure_function_1d(signal, tau_max) % 结构函数法计算一维信号分形维数 % signal : 表面轮廓, 行向量 % tau_max : 最大滞后距离, 一般取信号长度的 1/4 到 1/2 signal = signal(:)'; N = length(signal); if nargin < 2 tau_max = floor(N / 4); end taus = 1:tau_max; S = zeros(size(taus)); for idx = 1:length(taus) tau = taus(idx); diff_sq = (signal(1:end-tau) - signal(1+tau:end)).^2; S(idx) = mean(diff_sq); end % 对数域线性拟合 x = log(taus(:)); y = log(S(:)); p = polyfit(x, y, 1); H = p(1) / 2; D = 2 - H; end结构函数法的代码比功率谱法简单，但参数陷阱藏在两处：滞后范围和信号平稳性。tau_max取信号长度的 1/4 是常用经验值，因为 τ 接近信号长度时，参与平均的差分对只有个位数，方差急增，对数域散点会严重发散。另一个问题是信号如果有全局趋势（比如整个表面呈波浪形），结构函数在小 τ 区间反映的是局部分形特征，大 τ 区间反映的却是宏观形状，两者混合后拟合斜率会出现明显的拐点。解决办法是对信号做去趋势处理，或者只取小 τ 段（比如tau_max设为信号长度的 1/8）进行拟合。
二维表面的结构函数法和一维类似，但需要按方向计算。可以按行计算结构函数后取平均，也可以沿两个正交方向分别计算再取平均。对于各向同性表面，两种做法结果一致；对于各向异性表面，分方向算出的维数差异本身就是有意义的特征参数。
3.3 三种方法的数学关系和适用边界
三种方法并非独立。从理论上看，分形布朗运动的功率谱斜率、结构函数斜率、盒计数维数之间存在固定换算关系。对于一维分形信号：DBC 维数 D_dbc、功率谱斜率 β、结构函数斜率 α 之间满足D_dbc ≈ D_psf ≈ (5 - β)/2 ≈ (4 - α)/2。实际上，由于不同算法对尺度的加权方式不同，计算值会有系统性偏差。DBC 对表面局部奇异性更敏感，倾向于给出更高维数；功率谱法受频段截取影响大，频谱窗户选择不同偏差可达到 0.15；结构函数法对噪声最不敏感，但需要信号满足平稳增量条件。
什么时候用哪种方法？我的建议是：图像纹理分类首选 DBC，因为不需要标定物理尺度，直接算灰度即可；表面粗糙度分析且有标定高度数据时，结构函数法优先，因为它能给出 Hurst 指数，工程上可以用 H 值判断表面是偏“记忆性”（H > 0.5）还是偏“反记忆性”（H < 0.5）；需要分析周期性纹理或分形方向性时，功率谱法不可替代。
4. 参数调试与边界处理：三个必调参数和三条血泪经验
4.1 盒子尺寸序列：范围怎么定才不虚高
DBC 最常被问的参数就是box_sizes怎么给。粗糙的做法是直接[2 4 8 16 32 64]，但实际图像尺寸可能是 512、1024 甚至更大，这时小尺寸端的s=2会造成严重的高估。原因是：当盒子边长只有 2 个像素时，灰度方向被分成了G/s层，s=2意味着最多把 256 级灰度切成 128 层，而真实灰度动态范围远小于 256，导致大量空盒子被计入统计。
我处理这个问题的方法是：box_sizes从 4 或 6 开始，而不是从 2 开始，最大尺寸取图像短边长度的 1/4 到 1/8，并且在拟合时画散点图观察线性程度，只保留线性良好的区间。一个可靠的经验法则是让盒子尺寸序列长度大于等于 5 个点，如果只有 3 个点，拟合斜率受单点噪声影响极大。另外，盒子尺寸不必非用 2 的幂次，可以用round(linspace(4, max_size, 8))生成等差数列，但 2 的幂次在后续对比多尺度特征时更容易对齐。
提示：当盒子最大尺寸超过图像短边的 1/2 时，网格数量会小于 4 个，DBC 的统计根基就碎了。所有软件包里写死
[2,4,8,...,64]的默认值，在 256×256 图像上都会翻车。
4.2 灰度级裁剪上下限：图像动态范围必须铺满
DBC 的盒子高度与灰度级数 G 直接相关。代码里我取了整幅图的max_gray和min_gray做归一化，这个做法在图像整体对比度高时没问题，但遇到少量极端像素（比如金属断面的高亮点或孔洞阴影）时，整个动态范围被拉得过大，真正的纹理细节被压缩到很少的灰度层，导致盒计数在小尺度下几乎不随盒子尺寸变化，维数被显著低估。
更稳的做法是截取灰度分布的 1% 到 99% 分位点作为归一化范围，把极端像素都裁剪掉。mat2gray函数本身只做线性归一化，不能解决这个问题。我通常先跑prctile(I(:), [1 99])看一下灰度分布，如果 1 分位到 99 分位的跨度小于整幅图像动态范围的三分之一，就需要谨慎处理极端值，而不是直接归一化。
4.3 最小二乘拟合的截断策略：不能全段照单全收
三种方法最后都落在对数域线性拟合上，最常见的坑是：直接对全部数据点做polyfit，但实际对数域散点往往两端弯曲——小尺度端受像素离散化影响，大尺度端受样本有限性影响，中间一段才是线性区域。我处理这些数据时，会先用loglog画图，找到中间线性区间，再用polyfit对截取后的区间进行拟合。这个方法虽然看起来像“玄学”，但后来发现它确实能稳定地避免两端发散的问题。
一个量化判断标准是看拟合残差标准差：对残差值超过某个阈值的数据点直接剔除，然后重新拟合。Matlab 中可以用robustfit代替polyfit做鲁棒拟合，它能自动降低离群点权重，实测效果比手动截断更稳定。实际计算中，我一般是这样处理的：先用loglog画出曲线观察线性区间，或者直接使用robustfit避免手动截断。
功率谱法在处理窗函数时的边缘泄漏问题也很常见：图像边缘灰度对不齐会产生虚假高频分量，导致高频段功率谱抬升、斜率变缓。解决方法是先做镜像延拓或加 Hann 窗再做 FFT。镜像延拓的实现是I_padded = [I fliplr(I); flipud(I) fliplr(flipud(I))]，这样边缘不连续变成了对称过度，可以显著改善高频段斜率稳定性。
4.4 三个方法的翻车案例
A 同学做 SEM 图像形态分析，直接用 DBC 默认参数算出维数 2.6，明显偏高。检查发现图像里有一个标签残影，把图像的动态范围拉宽，真实纹理的灰度差被压扁。后来改用分位数裁剪，维数稳定在 2.1 附近，反而与手动估测一致。
另一个典型案例是对磨损表面做结构函数计算，tau_max 取信号长度的 1/2，导致大 τ 区间统计失真，拟合出现拐点。把 tau_max 改为 1/8，曲线恢复线性。实际案例中，磨损表面轮廓在特定方向上表现出多尺度特性，不同尺度的斜率变化本身验证了表面各向异性。这提醒我结构函数法当 tau 接近信号长度时需要谨慎解释，最好是加一个斜率滑窗检测。
功率谱法做表面粗糙度分析时，如果原始表面包含强周期性波纹，功率谱在对数坐标下会出现一个尖峰，整个拟合直线被这个尖峰拉偏。处理办法是先做空间滤波去除周期分量，或者拟合时手动跳过尖峰对应的频率区间。综合来看，唯一有效的方法是分频段拟合——把频谱分成低频段和高频段分别拟合，
本文还有配套的精品资源，点击获取
