简介：面向图像处理、信号分析与复杂系统研究，这套MATLAB分形维数计算资源包提供了差分盒维数、功率谱和结构函数三种主要算法的代码实现，能够帮助解决非规则几何对象难以量化建模的问题。压缩包内共5个m文件，整体仅3KB，代码紧凑地覆盖了不同尺度盒子的统计与log-log斜率拟合、基于快速傅立叶变换的频谱成分提取，以及多尺度结构函数变化特征的计算，可直接在MATLAB环境中运行并对照算法原理逐段阅读。资源已有2911人学习，适合需要复现分形实验或开展论文仿真的学生和研究者，既可将其作为快速计算工具，又能够借助简洁的函数拆分扩展出面向自身样本数据的分析流程。通过运行和调试这些脚本，读者能够深入理解分形维数的计算逻辑与自相似性量化方式，并迁移至图像纹理识别、表面粗糙度分析等实际科研任务，从而节省大量底层编码时间。 我叫不上“分形维数”这个词的，大多是做粗糙表面分析、信号特征提取或者图像纹理识别的人。可一旦开始搜Matlab实现，你会发现网上的代码版本五花八门，差分盒维数、功率谱、结构函数三种方法各说各话，跑出来的结果还经常对不上。我最初接触这个方向是为了做表面轮廓的粗糙度量化，被这些算法折腾了将近两周，踩了不少坑，也把每种方法的脾气摸了个大概。这篇东西就把我在实际计算中用到的Matlab实现、参数坑和选择逻辑一次性说清楚。
1. 为什么分形维数不是“一个数”，而是一组度量方式
先解决一个最根本的疑惑：为什么同一组数据，用差分盒维数算出来是2.31，用功率谱法算出来却成了1.65？这不是代码写错了，而是你根本没搞清楚自己该用哪种定义。
分形维数本质上量化的是“表面或曲线在不同尺度下表现出的复杂程度”。它不像长度、面积那样有一个绝对标准答案，不同的度量方式从不同角度刻画这种复杂度。你可以把分形维数理解成描述“细节填充程度”的指标——同样是100米长的海岸线，你用1公里尺子量是一个长度，用1米尺子量又是一个长度，尺子越细量出来的总长度越长。分形维数衡量的就是这种“长度随尺子缩小而增长的速度”。
对于二维图像或粗糙表面，常用的分形维数有三种定义路径：
- 盒维数（Box-counting / Differential Box-counting）：用不同大小的盒子去覆盖图像表面，统计非空盒子数随盒子尺寸的变化率。它衡量的是几何覆盖率，适合处理图像灰度表面这类二维数据。
- 功率谱法（Power Spectral Density method）：对数据做FFT变换后，统计功率谱密度与频率在双对数坐标下的斜率。它衡量的是能量在不同空间频率上的分布规律，本质上是自仿射分形特征的频域表达。
- 结构函数法（Structure Function method）：统计数据在不同滞后间隔下的差分二阶矩与滞后量的关系，从空间域的角度刻画自仿射分形的尺度行为。
这三种方法不是互相替代的关系，而是分别从覆盖率、频域能量、空间差分三个视角评估同一组数据的分形特性。实际应用中，表面轮廓分析更常用结构函数法，图像纹理分析更常用差分盒维数，而信号与噪声特征提取则首选功率谱法。
2. 差分盒维数的Matlab实现：图像粗糙度分析的主力算法
差分盒维数（Differential Box-Counting，DBC）是我最早接触也最常用的一种方法。它特别适合处理灰度图像的分形特征计算，很多材料科学中表面形貌分析、医学图像纹理分类的论文都基于这种方法。
2.1 算法原理与核心步骤
差分盒维数的思路是：把一张M×N的灰度图像想象成一个三维曲面，像素坐标(x,y)表示位置，灰度值z表示高度。然后用边长为s的盒子去覆盖这个三维曲面，统计需要多少个盒子才能完全覆盖曲面上的所有点。
具体步骤是这样的：
- 将图像划分为大小为s×s的网格块，假设图像尺寸为M×N，则每个网格块对应一个s×s的像素区域。
- 对于每个网格块，统计该区域内灰度值的最大值和最小值。灰度值范围为[0,255]，将其除以盒子边长s，得到灰度层数G = 255 / s。每个网格块内需要覆盖从最小灰度到最大灰度之间的盒子数n(i,j) = ceil(max/s) - ceil(min/s) + 1。
- 将所有网格块需要的盒子数累加，得到总盒子数N(s)。
- 改变盒子边长s的大小，重复上述过程，得到一系列(s, N(s))数据对。
- 对log(1/s)和log(N(s))做线性回归，斜率就是分形维数D。
分形维数的物理意义是：如果斜率越接近3，说明表面越“粗糙”，细节越丰富；越接近2，说明表面越“平滑”，纹理越简单。实际材料表面分形维数通常落在2到3之间。
2.2 直接可跑的Matlab代码
下面是我调试过多次、实际可用的差分盒维数函数：
function D = dbc_fractal_dimension(I) % I: 输入灰度图像矩阵，double类型，范围[0,255] % D: 计算得到的分形维数 I = double(I); [M, N] = size(I); if M ~= N error('图像必须是正方形'); end % 盒子尺寸序列，取2的幂次 L = floor(log2(min(M, N))); s = 2.^(1:L); % 从2到2^L num_scales = length(s); log_inv_s = zeros(num_scales, 1); log_Ns = zeros(num_scales, 1); for k = 1:num_scales box_size = s(k); % 将图像划分为网格 grid_size = M / box_size; % 对每个网格块计算盒子数 nr = floor(M / box_size); nc = floor(N / box_size); total_boxes = 0; for i = 0:nr-1 for j = 0:nc-1 % 提取当前网格块的灰度范围 block = I(i*box_size+1:(i+1)*box_size, j*box_size+1:(j+1)*box_size); min_val = min(block(:)); max_val = max(block(:)); % 计算该块需要的盒子数，注意+1的处理 n_ij = ceil(max_val/box_size) - ceil(min_val/box_size) + 1; total_boxes = total_boxes + n_ij; end end log_inv_s(k) = log(1/box_size); log_Ns(k) = log(total_boxes); end % 最小二乘线性拟合 p = polyfit(log_inv_s, log_Ns, 1); D = p(1); end2.3 实际计算时的关键细节与坑点
跑过这个算法的人多半遇到过拟合斜率不稳定、不同尺寸段结果差异大的问题。我总结一下最有价值的几个经验：
盒子尺寸序列的选择直接决定结果可靠性。我一开始直接用s从2到min(M,N)等间隔取，结果小尺寸段噪声特别大，大尺寸段又因为网格块数量太少导致统计失真。后来改成s取2的幂次（2,4,8,16...），拟合结果稳定多了。因为盒子尺寸翻倍对应着尺度翻倍，在双对数坐标下数据点分布更均匀，回归结果更可信。
小尺寸下线性拟合偏差大。盒子尺寸为2的时候，灰度分层只有几十层，统计误差非常大。建议在实际拟合时剔除s=2的数据点，从s=4开始拟合，斜率会更接近理论值。
图像尺寸必须是2的幂次倍数。如果不满足，边缘会有大量无法划分的网格块，导致盒子数统计偏低。最简单的做法是先用imresize把图像缩放到256×256或512×512。实测下来512×512的结果最稳定，256×256偏小约0.05左右。
灰度分层用ceil和floor混合处理时需要统一口径。有些实现直接用floor(max/s) - floor(min/s)，但这样会少算一个盒子，导致结果整体偏大0.1左右。用ceil(max/s) - ceil(min/s) + 1是公认比较准确的处理方式。
3. 功率谱法：一维信号分形维数计算的频域路径
功率谱法是从频域计算分形维数的方法，特别适合处理一维信号，比如表面轮廓曲线、时间序列、粗糙度测量数据。它的数学基础是：自仿射分形信号的功率谱密度满足幂律关系S(f) ∝ f^(-β)，其中β与分形维数D存在确定关系。
3.1 原理推导与公式使用
对于一维信号，功率谱密度函数和分形维数之间存在这样的关系：
β = 5 - 2D_h
其中D_h是赫斯特指数相关的分形维数。在标准的一维自仿射分形定义下，分形维数D与功率谱指数β的关系是：
D = (5 - β) / 2
也就是说，只要在双对数坐标下拟合功率谱的斜率β，就能直接算出分形维数。β的范围通常在1到3之间。β越接近3，D越接近1，信号越平滑；β越接近1，D越接近2，信号越粗糙。
功率谱法的核心优势是物理意义清晰——它直接告诉你能量如何随频率衰减。如果信号中存在周期性成分，功率谱上会出现明显的峰值，这会干扰斜率拟合，需要先做去趋势或滤波处理。
3.2 完整Matlab实现
function D = psd_fractal_dimension(x, fs) % x: 输入一维信号 % fs: 采样频率，可不传入，默认1 % D: 分形维数 if nargin < 2 fs = 1; end % 去除线性趋势 x = detrend(x(:), 'linear'); % 做FFT并计算功率谱密度 N = length(x); X = fft(x); Pxx = abs(X(1:floor(N/2)+1)).^2 / (fs * N); freq = (0:floor(N/2)) * fs / N; % 去掉直流分量 Pxx = Pxx(2:end); freq = freq(2:end); % 忽略低频段，防止趋势残留影响 mask = freq > 0.01 * fs; freq_fit = freq(mask); Pxx_fit = Pxx(mask); % 双对数回归 p = polyfit(log(freq_fit), log(Pxx_fit), 1); beta = -p(1); D = (5 - beta) / 2; end3.3 频域方法的典型问题与处理技巧
用功率谱法时最容易翻车的几个点，我一个个说：
FFT之前必须做去趋势。如果信号存在线性趋势，功率谱在低频段会严重偏离幂律关系，导致斜率估计偏低。detrend函数做线性去趋势是最基础的，如果还有高阶趋势，建议先做多项式拟合再去除。这个坑我踩过——有一组实验结果拟合出的分形维数只有1.1，怎么都不对，去掉趋势后变成1.72，完全合理。
频率段的选择范围要控制好。从直流分量往后开始拟合是非常糟糕的做法，因为低频段通常包含了趋势残留或仪器漂移。建议从特征频率的1/100处开始拟合，到奈奎斯特频率结束。同时，接近奈奎斯特频率的高频段也会因为FFT泄漏产生偏差，实际拟合时可以只取到fs/4甚至fs/8。
FFT点数不够时结果偏小。如果信号长度只有几百个点，功率谱在低频段的点太少，斜率拟合不稳定。建议信号长度至少2048点以上，数据不够可以用插值或重复测量补齐。
窗函数的影响需要评估。直接对原始信号做FFT会引入频谱泄漏。加汉宁窗可以抑制泄漏，但也会改变功率谱的斜率。实测下来加汉宁窗后β会略微偏大（分值维数偏小0.02-0.05）。如果只是做相对比较，可以忽略；如果需要绝对精度，建议不加窗或用平顶窗，并在论文中注明。
4. 结构函数法：粗糙表面分形分析的最稳健选择
结构函数法是分形维数计算中抗噪性能最好、对数据长度要求最低的方法之一。很多做表面形貌分析的老工程师更信任这种方法，因为它直接从空间域计算，避免了FFT固有的边界效应和频谱泄漏问题。
4.1 结构函数法的数学原理
结构函数法的核心是计算数据在不同滞后距离下的差分二阶矩。对于一条轮廓曲线z(x)，定义结构函数为：
S(τ) = E{|z(x+τ) - z(x)|²}
其中E表示数学期望，τ为滞后距离。对于自仿射分形轮廓，结构函数与滞后距离满足幂律关系：
S(τ) ∝ τ^(2H)
其中H为赫斯特指数，取值范围0到1。一维轮廓的分形维数与赫斯特指数的关系为：
D = 2 - H
因此，对log(τ)和log(S(τ))做线性回归，得到斜率2H，再算出D。
结构函数法对应的物理场景非常直接：滞后距离越小，差分值越小；滞后距离越大，差分值越大。表面越粗糙，不同滞后尺度上的差分值差异就越小（H越小，D越大）；表面越平滑，差分值随滞后距离增长越快（H越大，D越小）。
4.2 Matlab代码与参数设置
function D = structure_function_fd(z, max_tau_ratio) % z: 输入轮廓数据，一维数组 % max_tau_ratio: 最大滞后距离占数据长度的比例，默认0.25 % D: 分形维数 if nargin < 2 max_tau_ratio = 0.25; end z = z(:)'; N = length(z); max_tau = round(N * max_tau_ratio); tau_values = 1:max_tau; S_values = zeros(size(tau_values)); for tau = tau_values diffs = z(1+tau:end) - z(1:end-tau); S_values(tau) = mean(diffs.^2); end % 只取S>0的数据点做拟合 valid = S_values > 0; tau_fit = tau_values(valid); S_fit = S_values(valid); % 双对数线性回归 p = polyfit(log(tau_fit), log(S_fit), 1); H = p(1) / 2; D = 2 - H; end4.3 结构函数法的关键经验
最大滞后距离必须限制。如果τ接近数据长度N，参与平均的差分点太少，S(τ)的方差急剧增大，拟合结果完全失控。文献中建议最大滞后取数据长度的1/4到1/5，我在实际使用中一般取1/5，稳定性更好。如果数据长度只有500点，最大滞后取100，能用的拟合段大概是10到100，回归效果还可以。
数据采样密度决定可用尺度范围。结构函数法在小滞后段（τ小于采样间隔对应空间分辨率）没有物理意义。如果轮廓数据是等间距采样的，τ=1对应的空间尺度就是采样间隔。仪器噪声通常集中在最小的几个滞后尺度上，建议实际拟合时跳过τ=1到τ=3，从中间尺度开始。
三种方法可以用标准分形曲线验证。我当时用Weierstrass-Mandelbrot函数生成了已知分形维数的人造轮廓来验证程序正确性。例如生成D=1.5的W-M函数曲线，分别用结构函数法和功率谱法计算，结构函数法误差在0.02以内，功率谱法误差在0.05以内，差分盒维数需要二维化处理不太适合一维验证。
5. 三种方法的横向对比与选型建议
用对了方法，分形维数才有意义。这里我按实际经验把三种方法做了对比，并提供选型建议。
| 对比维度 | 差分盒维数 | 功率谱法 | 结构函数法 |
|---|---|---|---|
| 适用数据类型 | 二维灰度图像 | 一维信号/图像表面 | 一维轮廓/粗糙度曲线 |
| 计算原理 | 几何覆盖 | 频域幂律拟合 | 空间差分统计 |
| 抗噪性能 | 中等 | 较弱 | 较强 |
| 对数据长度要求 | 图像≥128×128 | 信号≥2048点 | 数据≥500点 |
| 对边界效应敏感度 | 中等 | 高 | 低 |
| 典型误差范围 | 0.03~0.08 | 0.05~0.12 | 0.01~0.05 |
| 计算速度 | 慢（双重循环） | 快 | 中等 |
从实际工程角度看，我给出这样几个选型建议：
- 做图像纹理分析，就用差分盒维数。虽然计算慢，但它对二维数据最友好，纹理分类效果也经过大量论文验证。
- 做一维信号分析且数据量充足，功率谱法和结构函数法都行。但如果你关心的是信号在不同频段的能量分布特征，功率谱法能提供额外信息，更适合。
- 做表面轮廓或粗糙度分析，首选结构函数法。抗噪能力强，对边界不敏感，数据要求也较低。
还有一个容易被忽略的经验：三种方法计算同一对象的理论分形维数结果并不完全一致，这很正常。论文中必须明确标注用的是哪种方法，否则审稿人第一件事就是质疑你。我在写报告时会同时算两种方法，如果结果在0.2以内一致，就说明分析可靠；如果差异超过0.3，就说明数据本身不符合自仿射分形的假设，需要检查数据来源。
6. 我踩过的几个坑与最后一组建议
最后分享几个实操层面的教训，这些在教科书和论文里基本找不到，但对实际项目帮助很大：
FFT的功率谱密度不同定义会导致结果差很多。Matlab的periodogram函数和直接做的fft算出来的Pxx绝对值差一个比例因子，但log-log拟合的斜率不受影响。真正有影响的是单边谱还是双边谱的选取——如果用双边谱，频率范围会扩大两倍，但功率谱密度值也会相应调整，斜率不变。所以这个不算坑，真正的问题是：如果你在低频段不剔除直流和趋势项，斜率会假性偏小（β偏大，D偏小0.2甚至更多）。
差分盒维数在图像分辨率变化时结果会漂移。同一个表面，用100×100像素的图和200×200像素的图计算，结果差异可能达到0.1。这不是算法问题，而是因为图像分辨率改变了可观察的最小尺度范围。做系列实验时一定要固定图像采集参数，要么所有图像都缩放到相同尺寸再比较。
数据处理流程比算法本身更影响结果。很多人在分形维数计算上纠结算法细节，却忽略了前置处理。去噪、去趋势、滤波、归一化这些步骤，每一项都会影响最终结果。我的建议是：在一组对比实验中，所有数据必须走完全相同的预处理流程，哪怕预处理参数稍微不合理——至少能保证相对比较的有效性。
回想起来，分形维数计算的难点从来不是算法有多复杂，而是你知不知道每种方法的前提假设是什么、数据的物理特征是什么、两者匹不匹配。代码跑起来只是第一步，理解数据背后的尺度行为才是真正能支撑你写出可靠分析结果的核心能力。
本文还有配套的精品资源，点击获取
