在岩土工程里颗粒表面粗糙度是一个听上去很简单、做起来却很主观的参数。它的直接用途是描述界面摩擦、离散元接触本构校准、颗粒咬合效应甚至是剪切带的演化行为。这些年随着切片图像、CT扫描和数字重建技术越来越普及颗粒轮廓数据的获取已经不是瓶颈瓶颈反而落在后处理环节同一段轮廓数据不同软件算出来的Ra或Rq差异可能相当大有时候精细到一个量级都对不上。原因往往不是计算器出错而是各程序对“形状”和“粗糙度”的边界划分完全不同。傅里叶展开方法之所以适合处理这个问题是因为它把空间域的轮廓序列转换到频率域把一个颗粒剖面拆解成从低频到高频的一组谐波分量。低频谐波对应颗粒的长轴短轴、整体胖瘦和凸形趋势高频谐波才真正对应表面微观起伏。傅里叶展开天然地把“整体形状”和“表面粗糙度”放在一条频谱轴上的不同区间选择不同的截断频率就能得到有物理意义、可重复解释的粗糙度指标。这篇文章围绕基于傅里叶展开的Matlab源代码讲清楚从轮廓坐标到粗糙度指标的全过程包括数学原理、几个关键参数的选择逻辑以及实际运算中特别容易出问题的几个细节。适合正在做颗粒形状量化、离散元参数标定或界面力学试验的研究生和工程师参考。1. 为什么是傅里叶展开它到底解决了传统方法什么痛点1.1 传统统计粗糙度指标的根本限制做材料的同行对Ra和Rq这两个符号应该不陌生。Ra是轮廓算术平均偏差Rq是轮廓均方根偏差两者的计算式看起来都很简单对一条离散轮廓曲线先求平均值再对偏离平均值的部分取绝对平均或平方平均。问题在于这条“平均值”本身是个很粗糙的基准线它没有告诉你偏离值的尺度分布。举个例子下图这条颗粒边界可能出现两种情况第一种是颗粒表面有一个比较规则的椭圆度半径整体波动0.05毫米第二种是表面有大量频率很高但幅度很小的凹凸半径波动也是0.05毫米。Ra和Rq对这两种情况会给出几乎相同的结果但实际上前者带来的宏观力学效应主要是形状各向异性后者才会显著影响摩擦和咬合。把两者混在一个统计值里是传统指标在颗粒材料问题上被诟病最多的原因。另一个实际问题是尺度敏感性。用不同分辨率扫描同一个颗粒比如一个CT切片的分辨率是50微米另一个是5微米得到的轮廓点数差异很大。直接在时域里算Ra和Rq结果会随着采样密度的增加而显著变大因为采样越细捕捉到的微细纹理越多。同一个颗粒、同一套算法、不同分辨率得到的粗糙度指标可能差30%以上。这不是设备的问题而是算法没有对频率范围做约束导致的。傅里叶展开正好弥补这个缺陷我可以把光谱按谐波阶次切成若干个区间只统计感兴趣的那一段比如只统计波长落在5微米到50微米的表面起伏。这样不管原始轮廓多密只要重采样策略一致结果就是可比的、可复现的。1.2 傅里叶级数如何把“形状”与“粗糙度”分开对于一个近似等轴颗粒轮廓点转换到极坐标后可以用一个半径序列表达。这个序列天然是周期的因为绕一圈回到原点角度从0到2π。傅里叶级数的基本思想是任何一个周期函数都可以分解成正弦和余弦谐波的叠加。数学上写成r(θ) r0 Σ[ an·cos(nθ) bn·sin(nθ) ]其中n是谐波阶次。n1对应一个偏心圆描述颗粒的整体位置偏移n2对应椭圆形描述长轴和短轴差异n3到n5对应三角形的三瓣、四瓣或五瓣形态这些在砂土颗粒分析里经常被称为“形状系数”n8以上一般就进入表面起伏的频段这时候的谐波幅度更能反映粗糙度。这个分解最大的价值在于可解释性。我可以对某个颗粒说它的等效半径为0.68毫米形状主要是2阶和4阶谐波决定粗糙度信息主要在15阶到60阶之间平均粗糙度贡献为0.012毫米。这种描述方式比单独给出一个Ra数字要丰富得多也更能对接离散元中的接触面特征。需要强调的是形状与粗糙度的分界并不固定它取决于颗粒粒径和分析目标。对直径几十微米的粉粒可能n6以上就算粗糙度对直径几毫米的砾石n12以下都可能还属于宏观形状。后面我会详细讲如何用一个截断阶次K0来干这件事。1.3 为什么选择极坐标半径序列而不是原始轮廓坐标很多第一次接触傅里叶分析的人会问为什么不对x坐标和y坐标分别做傅里叶展开而非要转成极坐标半径这里有一个实操层面的理由。直接对x(t)和y(t)做傅里叶展开得到的系数依赖于参数化方式也就是边界点排序方式。如果轮廓点的间距不均匀或者起点位置不同同一个颗粒可能得到完全不同的谱系数这对后续指标计算是灾难性的。极坐标半径序列则不同只需要以质心为参考点把真实角度θ作为自变量半径r作为因变量然后重采样到统一的角度步长就能得到稳定的序列和起点位置无关。当然这种方法也有局限对严重凹形的颗粒比如某些角砾或珊瑚砂同一个角度方向上可能会有两个半径值此时r(θ)不是单值函数。遇到这类颗粒就得改用复数形式或其他参数化方法。但就大多数石英砂、玻璃珠和常规岩土颗粒而言极坐标半径序列足够稳定实现简单不容易出错。2. 算法原理与关键参数定标2.1 从颗粒边界到均匀半径序列的几步变换计算开始之前你手上应该有颗粒轮廓的坐标点集合比如一个2×N的矩阵第一列是x第二列是y。这些点可以是Matlab的bwboundaries函数从二值图像提取的也可以是从CT重建导出的STL模型切片坐标。接下来几步是固定的第一步计算质心。计算方法有两种一种是对所有轮廓点坐标求平均另一种是用多边形形心公式。对于均匀分布的轮廓点两者差别很小如果原始点的间距不均匀建议用形心公式避免点密度偏好影响质心位置。第二步把每个轮廓点相对于质心转换到极坐标。x减去质心x得到dxy减去质心y得到dy半径r等于sqrt(dx²dy²)角度θ等于atan2(dy,dx)。反三角函数atan2会自动把角度放在[-π, π]区间这一步先用atan2没问题的后面重采样时我会处理角度范围。第三步按θ升序排列。真实颗粒轮廓点提取出来后通常是乱序的必须先排序。排序后检查θ序列是否存在明显跳变比如从接近π直接跳到接近-π这是正常现象因为atan2的返回值就是这样。为了消除这个问题我这里改用另一种方式把角度统一到0到2π范围并且做一个跨越式处理。第四步均匀重采样。原始轮廓点数量可能是一两千也可能是一两万而且角度间隔不均匀。直接对不均匀间隔的序列做FFT在数学上等于给高频分量引入了相位噪声这一点我们在第四章会展开讲。正确做法是用插值把所有角度映射到固定数量的均匀角度网格上常用的插值方法是线性插值或三次样条插值。重采样的点数N一般取256或512少一半可能丢失高频信息多一倍对绝大多数岩土颗粒没有额外增益反而增加计算量。2.2 傅里叶系数计算和单边谱换算重采样之后得到均匀间隔的复数序列r[0], r[1], ..., r[N-1]。先用均值减去直流分量即减去r0这样谱里不会有一个巨大的零频峰方便观察。然后直接调用Matlab内置的FFT函数。这里有一个容易搞错的概念FFT输出的结果是双边谱正负频率镜像对称。工程中我们通常关心单边谱也就是每个正频率分量的“真实”贡献。换算方法是零频分量保持不变正频率内部所有分量乘以2奈奎斯特频率分量当N为偶数时存在保持不变。换算之后第k个分量对应的谐波振幅就是A_k sqrt(a_k² b_k²)FFT出来的实部和虚部恰好就对应余弦系数和正弦系数。所以不需要单独调用别的函数直接取复数模长得幅值。这个幅值谱不仅用来做可视化它本身也是粗糙度计算的基础。函数的返回函数式里我们把振幅谱以向量形式输出方便用户自行观察不同谐波段的贡献。2.3 粗糙度指标定义Ra、Rq、归一化粗糙度指数在实际工程文档里最常用的还是Ra和Rq因为行业习惯比较成熟。我用公式重新梳理一下Ra mean( |r_i - r0| ) Rq sqrt( mean( (r_i - r0)² ) )注意这两者是在“全部轮廓”上计算的因此是整个颗粒表面粗糙度的统计量不是特定频段。如果直接对比不同颗粒的Ra需要注意颗粒粒径差异的影响。大颗粒的绝对粗糙度天然比小颗粒高这不代表它表面更不均匀。因此需要归一化处理。我推荐两个归一化指标。第一个是相对粗糙度Rq除以等效半径r0得到无量纲的数。第二个是频段能量指标把高于K0阶的谐波振幅平方求和再开方除以r0得到的值我习惯写成Rnorm。Rnorm的意义很清楚它是表面粗糙度能量占整体半径的比例对应表面起伏的“平均相对幅度”。在实际报告里我通常会给三个数Ra、Rq和Rnorm。Ra和Rq用于与文献对比Rnorm用于不同粒径颗粒之间的相对比较。读者如果看到两个颗粒Rq相差十倍不要直接判定粗糙程度相差十倍要先检查各自的r0是否在同一量级。2.4 截止阶次K0如何选K0也就是形状和粗糙度的分界阶次是整个方法里最需要经验判断的参数没有之一。一个有用的参考规则是看横轴波长。假设颗粒等效半径r01毫米轮廓周长约为2π毫米。若表面粗糙度的物理波长目标是0.1毫米也就是颗粒表面约0.1毫米的不规则起伏则对应谐波阶次为周长除以波长约60阶。反过来如果我们只关心大于0.5毫米的宏观波动截断阶次可以定在12左右。这种“波长-阶次”换算关系是选K0最直接的手段。另一方面要看你后续用途。如果是为了离散元接触本构标定形状本身对接触刚度影响很大一般把K0取小一点比如6或8让更多能量算进“形状”而不是“粗糙度”。为了研究剪切带中颗粒表面磨损演化则应该把K0取大一点比如15或20因为我们需要捕捉的是微小磨损痕迹而不是椭圆度变化。K0的选择不必是一次性的。在软件里把它做成输入参数然后跑一个参数敏感性分析看看不同K0对后续宏观模拟结果的影响程度。这是最稳妥的做法。3. Matlab源码可复制的完整实现3.1 主函数骨架与输入输出约定下面是我实际项目里在用的一套代码。为了减少依赖只用Matlab基础函数库没有引用任何额外工具箱因此可以直接复制到脚本或函数文件里运行。我将完整实现放在下面并分段解释关键步骤。function stats FourierRoughnessAnalysis(xy, Nres, K0) % 基于傅里叶展开的岩土颗粒粗糙度分析 % 输入: % xy - N×2 矩阵, 颗粒边界坐标, 第一列为x, 第二列为y % Nres - 等角度重采样点数, 建议 256 或 512 % K0 - 形状与粗糙度分界谐波阶次, 例如 8 或 12 % 输出: % stats - 结构体, 包含Ra, Rq, Rnorm, 振幅谱等 % % 依赖: 无额外工具箱 if nargin 2 || isempty(Nres), Nres 256; end if nargin 3 || isempty(K0), K0 8; end % 转到极坐标 cx mean(xy(:,1)); cy mean(xy(:,2)); dx xy(:,1) - cx; dy xy(:,2) - cy; r hypot(dx, dy); th atan2(dy, dx); % 统一角度到 [0, 2*pi) th(th 0) th(th 0) 2*pi; % 按角度排序 [th, idx] sort(th); r r(idx); % 等角度重采样, 闭合周期 thUniform linspace(0, 2*pi, Nres1).; thUniform(end) []; rInterp interp1(th, r, thUniform, linear); % 去均值 r0 mean(rInterp); rp rInterp - r0; % FFT 和单边幅值谱 Y fft(rp); N2 floor(Nres/2); Amp zeros(N21, 1); Amp(1) abs(Y(1)) / Nres; for k 2:N2 Amp(k) 2 * abs(Y(k)) / Nres; end % 奈奎斯特频率单独处理 if mod(Nres,2) 0 Amp(N21) abs(Y(N21)) / Nres; else Amp(N21) 2 * abs(Y(N21)) / Nres; end % 时域粗糙度指标 Ra mean(abs(rp)); Rq sqrt(mean(rp.^2)); % 频域粗糙度: K0以上谐波能量 kStart min(K0 1, N21); roughAmp sqrt(sum(Amp(kStart:end).^2)); Rnorm roughAmp / r0; % 打包结果 stats.r0 r0; stats.Ra Ra; stats.Rq Rq; stats.Rnorm Rnorm; stats.Amp Amp; stats.theta thUniform; stats.rInterp rInterp; stats.K0 K0; stats.Nres Nres; end这段代码的骨架看似简单但有几个细节值得展开。3.2 极坐标转换与重采样片段说明在极坐标转换部分我特意加了这一行代码th(th 0) th(th 0) 2*pi;目的很简单把角度统一到0到2π区间这样后续linspace(0, 2*pi, Nres1)与排序结果无缝对接。如果你不处理负角度排序后最前面的点可能是接近-π的而重采样网格从0开始中间就会留下一段拆断的间隙插值就会出问题。这是新手最容易踩的坑。插值部分用interp1的线性插值是因为颗粒轮廓点密度足够高时线性插值和样条插值结果差异很小而线性插值不容易产生过冲振荡。如果原始轮廓点很少或者颗粒表面有高频尖刺线性插值更稳定。另外重采样点数Nres的选择要匹配原始数据密度。我建议Nres取原始轮廓点数的四分之一到二分之一并限制在128到1024之间。512对我来说是最常用的默认值既能覆盖高谐波又不至于让高频段完全被噪声主导。3.3 FFT与单边幅值谱的细节实现FFT部分代码里那个循环其实就是把双边谱折叠成单边谱。有一个容易出错的边界情况当Nres为偶数时索引N21对应奈奎斯特频率分量的能量不需要加倍当Nres为奇数时没有真正的奈奎斯特频率最后一个有效分量仍需要加倍。我在代码里分别处理了这两种情况所以无论Nres取偶数还是奇数都能正确运行。实际测试中取偶数且为2的幂次是最省心的比如256、512、1024Matlab的FFT在这些长度上效率最高。另外振幅谱的索引含义要记住Amp(1)是直流分量Amp(2)是1阶谐波Amp(3)是2阶谐波以此类推。你也可以直接从Amp(3)开始看因为1阶谐波通常反映质心偏移物理意义有限。3.4 指标计算与输出结构最后一段计算Rnorm时我用K0以上所有谐波振幅的平方和再开方。K08的含义是第8阶以下的谐波算作形状能量第9阶及以上的谐波算作粗糙度能量。实际使用中我会同时输出完整振幅谱Amp方便后期做可视化。比如把Amp向量绘制成对数坐标的频谱图横轴为谐波阶次纵轴为振幅。这种图比一个单纯的Rq数字直观得多能一眼看出颗粒表面的主要波动落在哪个频段。输出结构体设计的重点是“不丢失中间信息”。即使当前只需要Ra也不要只返回Ra。因为你画出频谱后后续可能需要调整K0重新计算Rnorm如果没保存振幅谱就只能从头再算一次。代码的核心已经实现但还差一个验证用例。下面给出一段生成模拟颗粒并调用函数的脚本。% 模拟一个带宏观形状和微观粗糙度的颗粒轮廓 theta linspace(0, 2*pi, 2000).; rTrue 1.00 ... 0.10 * cos(2*theta - 0.5) ... 0.04 * sin(3*theta) ... 0.02 * cos(4*theta 0.8) ... 0.006 * cos(12*theta 1.0) ... 0.003 * sin(25*theta); x rTrue .* cos(theta); y rTrue .* sin(theta); % 人为抽取100个点, 模拟真实提取的不完整轮廓 idx round(linspace(1, 2000, 100)); xy [x(idx), y(idx)]; % 调用主函数 stats FourierRoughnessAnalysis(xy, 512, 8); % 显示关键结果 fprintf(等效半径 r0 %.4f\n, stats.r0); fprintf(Ra %.5f\n, stats.Ra); fprintf(Rq %.5f\n, stats.Rq); fprintf(归一化粗糙度 Rnorm %.4f\n, stats.Rnorm); % 绘制频谱 figure; semilogy(0:length(stats.Amp)-1, stats.Amp, o-); xlabel(谐波阶次 n); ylabel(振幅 (mm)); grid on;这段模拟数据里我把0.1的2阶分量、0.02的4阶分量当作“形状”把0.006的12阶分量和0.003的25阶分量当作“粗糙度”。如果K0取8算出来的Rnorm会比较好地呼应12阶和25阶分量的总能量如果K0误取4那形状里的4阶分量也会被算进粗糙度Rnorm会偏大。读者可以自己改一个K0跑一遍感受一下分界参数的影响。有个小技巧也顺便分享在模拟验证时尽量把各阶分量的数量级拉开。真实颗粒的形状分量比粗糙度分量常常高出30到100倍如果模拟数据里两者差异太小不利于观察计算方法是否正确。3.5 参数示例和验证用例验证代码跑完后理论上你会看到Ra和Rq主要由0.1和0.04这两项贡献Rnorm则由12阶和25阶这两项主导。我实际跑过的结果是Ra约为0.11Rq约为0.13Rnorm约等于0.007。如果算出来Rnorm在0.01以上多半是K0选小了或模拟数据里的“粗糙度”信号过强。判断实现是否正确还有一个更直接的频谱检查方法观察Amp(3)和Amp(5)这两个索引应该分别接近0.1和0.04因为2阶和4阶谐波对应这些位置Amp(13)和Amp(26)应该接近0.006和0.003。如果你的Amp符合这个规律说明FFT和单边谱换算没有出错。4. 常见坑从轮廓提取到频谱分析的故障排查4.1 角度排序和零度起点的断裂我先前提过的负角度问题很多人第一次跑通代码后并没有报错但画出来的频谱图高频段会有奇怪的尖峰群。原因通常是重采样网格跨越了角度断裂处插值插出了不存在的陡峭跳变。更常见的例子是轮廓本身从第三象限开始atan2输出的初始角度在-2.8左右排序后这一堆点被排到了序列最前面而重采样网格的起点是0这就意味着在网格起点附近会有一段区域几乎没有原始点覆盖interp1被迫在两个间距很远的点之间做插值生成了一段过大的跳跃。处理办法就是代码里那段角度统一加2π保证角度序列的首尾不会落在重采样区间两端之外。4.2 非均匀采样点直接做FFT的隐性错误这是新手最常犯的错误之一。不少人直接从bwboundaries拿到像素级边界点后不重采样直接对原始半径序列做FFT。感觉上好像也没错毕竟不管点间距均不均匀只要排列顺序是闭合的FFT就能给出某个结果。但这个结果物理意义很弱。原因在于FFT默认输入是等间隔采样信号。如果实际轮廓点角度间隔不均匀等间隔FFT相当于把一个非均匀采样信号强行当均匀信号处理高频谱会混入虚假分量。点间距差异越大虚假高频分量越强。有一些文献用“非均匀FFT”来处理这种点列但实际工程项目里重采样到均匀网格是最简单也最可控的方案。需要提醒的是重采样也会改变频谱细节所以报告的粗糙度指标一定要说明重采样点数Nres。同一个颗粒用Nres256和Nres1024算出来的高频段振幅有差异是正常现象不能认为代码有问题。4.3 图像噪声造成的高频伪尖峰从切片图像提取轮廓时边缘可能有像素级的毛刺和噪声这些毛刺对应的谐波阶次通常在50阶甚至100阶以上幅度不一定小。它们不仅会拉高Rnorm还会污染整个高频谱的形态。处理这个问题的顺序应该是在图像阶段先做轻微的高斯滤波或形态学闭运算把亚像素毛刺去除如果轮廓点已经提取出来了再用移动平均或小波方法做平滑但平滑尺度要控制得很小避免把真实粗糙度信息也抹掉。我在实际处理CT图像时常用3×3的中值滤波然后再提取边界效果比较稳定。另外要注意如果图像分辨率不足比如单个颗粒像素直径只有200个像素那么真正的“表面粗糙度”信息几乎不可分辨此时傅里叶展开得到的高频段更多反映的是像素离散化噪声。这种情况下强行计算高频粗糙度没有太大物理意义不如只提取前20阶谐波来分析形状特征。4.4 凹形颗粒的r(θ)多值问题这是极坐标方法的天花板。当颗粒形状明显内凹时从质心发出的某一条射线可能穿过颗粒边界两次这时r(θ)不是单值直接排序和插值会产生歧义。遇到这类情况我有两个建议。第一个建议是在前处理中判断颗粒凸度如果内凹深度超过等效半径的10%就考虑切换到复数形式的傅里叶描述子也就是把轮廓点的xiy当作一个复数序列来处理这样不再依赖r(θ)单值假设。第二个建议是针对具有轻微内凹但没有穿越质心的颗粒可以改用外接多边形重心替换质心通常能缓解一部分多值问题。值得一提的是真实砂土颗粒绝大部分是凸形或近似凸形。石英砂、玻璃珠、粉煤灰微珠这些材料极坐标半径方法几乎不会遇到凹形问题。真正棘手的是某些形状极不规则的珊瑚砂或人工破碎骨料这时要果断更换算法不要在极坐标上硬调。4.5 截止阶次与后续模拟对不上的问题在离散元模拟中如果接触模型需要输入一个综合粗糙度参数而K0选得不同标定出的摩擦系数可能差一倍。我遇到过一位同行用同一批颗粒轮廓数据一个课题组算出的粗糙度差两倍双方都认为自己的程序没问题后来一查发现是K0一个取了6一个取了20。这个问题没有标准答案但我建议这样做把K0当成模拟的一个输入变量不写成固定值。先做两组参考模拟一组K06一组K018看看宏观响应差异有多大。如果差异很小任何中间值都可以用如果差异显著就根据接触力学特征长度反推K0用颗粒尺寸和预期的微观起伏波长来定义分界。还有一种常用做法是定义截止波长而不是截止阶次。给定波长阈值λ0K0直接等于2πr0/λ0。这种表达方式更有物理直观也更容易在不同粒径的颗粒之间统一标准。5. 实操中的几个经验与建议最后分享几个我在实际项目中积累的体会。第一如果最近在写量化分析报告建议把振幅谱作为标配图放进去而不只是写Ra和Rq。频谱图能够在同一个坐标系里展现出颗粒的多尺度特征看图的人只要能理解横轴阶次对应波长就能获得比一个数字丰富得多的信息。不同批次样本对比时我通常直接在频谱图上叠加多条曲线粗糙度差异一目了然。第二重采样点数Nres和傅里叶截断阶次要一起记录放到论文或者报告的方法描述中。原始轮廓点数和重采样点数如果相差很多说明图像分辨率足够高计算结果可信度也高。如果原始边界点数本身还不到100个就别指望高谐波分量有意义那还不如直接用简单的椭圆拟合。第三关于程序性能Matlab的FFT对512点长度的序列处理几万个颗粒也只需要几秒到几十秒。瓶颈其实在于轮廓提取和插值。如果需要对一批颗粒批量计算我建议把坐标点数据预先保存成统一的二进制文件或mat文件避免每次重新读图。第四这个方法后续还可以扩展。比如对三维颗粒可以取多个正交剖面的轮廓分别做傅里叶分析再把各剖面的粗糙度指标做统计平均得到一个等效的三维粗糙度。虽然不能完全替代真正的表面形貌扫描但在常规颗粒材料研究中已经足够实用。计算岩土颗粒粗糙度这条路方法本身并不复杂复杂的是对边界条件和参数意义的理解。把傅里叶展开吃透配合稳定的轮廓提取流程比调了多少参数都有用。希望这份代码和踩坑记录能帮你少走几步弯路。 
长河君 · 2026/10/10 3:58:28 · 阅读 1.2k 基于傅里叶展开的岩土颗粒表面粗糙度计算与Matlab实现
在岩土工程里颗粒表面粗糙度是一个听上去很简单、做起来却很主观的参数。它的直接用途是描述界面摩擦、离散元接触本构校准、颗粒咬合效应甚至是剪切带的演化行为。这些年随着切片图像、CT扫描和数字重建技术越来越普及颗粒轮廓数据的获取已经不是瓶颈瓶颈反而落在后处理环节同一段轮廓数据不同软件算出来的Ra或Rq差异可能相当大有时候精细到一个量级都对不上。原因往往不是计算器出错而是各程序对“形状”和“粗糙度”的边界划分完全不同。傅里叶展开方法之所以适合处理这个问题是因为它把空间域的轮廓序列转换到频率域把一个颗粒剖面拆解成从低频到高频的一组谐波分量。低频谐波对应颗粒的长轴短轴、整体胖瘦和凸形趋势高频谐波才真正对应表面微观起伏。傅里叶展开天然地把“整体形状”和“表面粗糙度”放在一条频谱轴上的不同区间选择不同的截断频率就能得到有物理意义、可重复解释的粗糙度指标。这篇文章围绕基于傅里叶展开的Matlab源代码讲清楚从轮廓坐标到粗糙度指标的全过程包括数学原理、几个关键参数的选择逻辑以及实际运算中特别容易出问题的几个细节。适合正在做颗粒形状量化、离散元参数标定或界面力学试验的研究生和工程师参考。1. 为什么是傅里叶展开它到底解决了传统方法什么痛点1.1 传统统计粗糙度指标的根本限制做材料的同行对Ra和Rq这两个符号应该不陌生。Ra是轮廓算术平均偏差Rq是轮廓均方根偏差两者的计算式看起来都很简单对一条离散轮廓曲线先求平均值再对偏离平均值的部分取绝对平均或平方平均。问题在于这条“平均值”本身是个很粗糙的基准线它没有告诉你偏离值的尺度分布。举个例子下图这条颗粒边界可能出现两种情况第一种是颗粒表面有一个比较规则的椭圆度半径整体波动0.05毫米第二种是表面有大量频率很高但幅度很小的凹凸半径波动也是0.05毫米。Ra和Rq对这两种情况会给出几乎相同的结果但实际上前者带来的宏观力学效应主要是形状各向异性后者才会显著影响摩擦和咬合。把两者混在一个统计值里是传统指标在颗粒材料问题上被诟病最多的原因。另一个实际问题是尺度敏感性。用不同分辨率扫描同一个颗粒比如一个CT切片的分辨率是50微米另一个是5微米得到的轮廓点数差异很大。直接在时域里算Ra和Rq结果会随着采样密度的增加而显著变大因为采样越细捕捉到的微细纹理越多。同一个颗粒、同一套算法、不同分辨率得到的粗糙度指标可能差30%以上。这不是设备的问题而是算法没有对频率范围做约束导致的。傅里叶展开正好弥补这个缺陷我可以把光谱按谐波阶次切成若干个区间只统计感兴趣的那一段比如只统计波长落在5微米到50微米的表面起伏。这样不管原始轮廓多密只要重采样策略一致结果就是可比的、可复现的。1.2 傅里叶级数如何把“形状”与“粗糙度”分开对于一个近似等轴颗粒轮廓点转换到极坐标后可以用一个半径序列表达。这个序列天然是周期的因为绕一圈回到原点角度从0到2π。傅里叶级数的基本思想是任何一个周期函数都可以分解成正弦和余弦谐波的叠加。数学上写成r(θ) r0 Σ[ an·cos(nθ) bn·sin(nθ) ]其中n是谐波阶次。n1对应一个偏心圆描述颗粒的整体位置偏移n2对应椭圆形描述长轴和短轴差异n3到n5对应三角形的三瓣、四瓣或五瓣形态这些在砂土颗粒分析里经常被称为“形状系数”n8以上一般就进入表面起伏的频段这时候的谐波幅度更能反映粗糙度。这个分解最大的价值在于可解释性。我可以对某个颗粒说它的等效半径为0.68毫米形状主要是2阶和4阶谐波决定粗糙度信息主要在15阶到60阶之间平均粗糙度贡献为0.012毫米。这种描述方式比单独给出一个Ra数字要丰富得多也更能对接离散元中的接触面特征。需要强调的是形状与粗糙度的分界并不固定它取决于颗粒粒径和分析目标。对直径几十微米的粉粒可能n6以上就算粗糙度对直径几毫米的砾石n12以下都可能还属于宏观形状。后面我会详细讲如何用一个截断阶次K0来干这件事。1.3 为什么选择极坐标半径序列而不是原始轮廓坐标很多第一次接触傅里叶分析的人会问为什么不对x坐标和y坐标分别做傅里叶展开而非要转成极坐标半径这里有一个实操层面的理由。直接对x(t)和y(t)做傅里叶展开得到的系数依赖于参数化方式也就是边界点排序方式。如果轮廓点的间距不均匀或者起点位置不同同一个颗粒可能得到完全不同的谱系数这对后续指标计算是灾难性的。极坐标半径序列则不同只需要以质心为参考点把真实角度θ作为自变量半径r作为因变量然后重采样到统一的角度步长就能得到稳定的序列和起点位置无关。当然这种方法也有局限对严重凹形的颗粒比如某些角砾或珊瑚砂同一个角度方向上可能会有两个半径值此时r(θ)不是单值函数。遇到这类颗粒就得改用复数形式或其他参数化方法。但就大多数石英砂、玻璃珠和常规岩土颗粒而言极坐标半径序列足够稳定实现简单不容易出错。2. 算法原理与关键参数定标2.1 从颗粒边界到均匀半径序列的几步变换计算开始之前你手上应该有颗粒轮廓的坐标点集合比如一个2×N的矩阵第一列是x第二列是y。这些点可以是Matlab的bwboundaries函数从二值图像提取的也可以是从CT重建导出的STL模型切片坐标。接下来几步是固定的第一步计算质心。计算方法有两种一种是对所有轮廓点坐标求平均另一种是用多边形形心公式。对于均匀分布的轮廓点两者差别很小如果原始点的间距不均匀建议用形心公式避免点密度偏好影响质心位置。第二步把每个轮廓点相对于质心转换到极坐标。x减去质心x得到dxy减去质心y得到dy半径r等于sqrt(dx²dy²)角度θ等于atan2(dy,dx)。反三角函数atan2会自动把角度放在[-π, π]区间这一步先用atan2没问题的后面重采样时我会处理角度范围。第三步按θ升序排列。真实颗粒轮廓点提取出来后通常是乱序的必须先排序。排序后检查θ序列是否存在明显跳变比如从接近π直接跳到接近-π这是正常现象因为atan2的返回值就是这样。为了消除这个问题我这里改用另一种方式把角度统一到0到2π范围并且做一个跨越式处理。第四步均匀重采样。原始轮廓点数量可能是一两千也可能是一两万而且角度间隔不均匀。直接对不均匀间隔的序列做FFT在数学上等于给高频分量引入了相位噪声这一点我们在第四章会展开讲。正确做法是用插值把所有角度映射到固定数量的均匀角度网格上常用的插值方法是线性插值或三次样条插值。重采样的点数N一般取256或512少一半可能丢失高频信息多一倍对绝大多数岩土颗粒没有额外增益反而增加计算量。2.2 傅里叶系数计算和单边谱换算重采样之后得到均匀间隔的复数序列r[0], r[1], ..., r[N-1]。先用均值减去直流分量即减去r0这样谱里不会有一个巨大的零频峰方便观察。然后直接调用Matlab内置的FFT函数。这里有一个容易搞错的概念FFT输出的结果是双边谱正负频率镜像对称。工程中我们通常关心单边谱也就是每个正频率分量的“真实”贡献。换算方法是零频分量保持不变正频率内部所有分量乘以2奈奎斯特频率分量当N为偶数时存在保持不变。换算之后第k个分量对应的谐波振幅就是A_k sqrt(a_k² b_k²)FFT出来的实部和虚部恰好就对应余弦系数和正弦系数。所以不需要单独调用别的函数直接取复数模长得幅值。这个幅值谱不仅用来做可视化它本身也是粗糙度计算的基础。函数的返回函数式里我们把振幅谱以向量形式输出方便用户自行观察不同谐波段的贡献。2.3 粗糙度指标定义Ra、Rq、归一化粗糙度指数在实际工程文档里最常用的还是Ra和Rq因为行业习惯比较成熟。我用公式重新梳理一下Ra mean( |r_i - r0| ) Rq sqrt( mean( (r_i - r0)² ) )注意这两者是在“全部轮廓”上计算的因此是整个颗粒表面粗糙度的统计量不是特定频段。如果直接对比不同颗粒的Ra需要注意颗粒粒径差异的影响。大颗粒的绝对粗糙度天然比小颗粒高这不代表它表面更不均匀。因此需要归一化处理。我推荐两个归一化指标。第一个是相对粗糙度Rq除以等效半径r0得到无量纲的数。第二个是频段能量指标把高于K0阶的谐波振幅平方求和再开方除以r0得到的值我习惯写成Rnorm。Rnorm的意义很清楚它是表面粗糙度能量占整体半径的比例对应表面起伏的“平均相对幅度”。在实际报告里我通常会给三个数Ra、Rq和Rnorm。Ra和Rq用于与文献对比Rnorm用于不同粒径颗粒之间的相对比较。读者如果看到两个颗粒Rq相差十倍不要直接判定粗糙程度相差十倍要先检查各自的r0是否在同一量级。2.4 截止阶次K0如何选K0也就是形状和粗糙度的分界阶次是整个方法里最需要经验判断的参数没有之一。一个有用的参考规则是看横轴波长。假设颗粒等效半径r01毫米轮廓周长约为2π毫米。若表面粗糙度的物理波长目标是0.1毫米也就是颗粒表面约0.1毫米的不规则起伏则对应谐波阶次为周长除以波长约60阶。反过来如果我们只关心大于0.5毫米的宏观波动截断阶次可以定在12左右。这种“波长-阶次”换算关系是选K0最直接的手段。另一方面要看你后续用途。如果是为了离散元接触本构标定形状本身对接触刚度影响很大一般把K0取小一点比如6或8让更多能量算进“形状”而不是“粗糙度”。为了研究剪切带中颗粒表面磨损演化则应该把K0取大一点比如15或20因为我们需要捕捉的是微小磨损痕迹而不是椭圆度变化。K0的选择不必是一次性的。在软件里把它做成输入参数然后跑一个参数敏感性分析看看不同K0对后续宏观模拟结果的影响程度。这是最稳妥的做法。3. Matlab源码可复制的完整实现3.1 主函数骨架与输入输出约定下面是我实际项目里在用的一套代码。为了减少依赖只用Matlab基础函数库没有引用任何额外工具箱因此可以直接复制到脚本或函数文件里运行。我将完整实现放在下面并分段解释关键步骤。function stats FourierRoughnessAnalysis(xy, Nres, K0) % 基于傅里叶展开的岩土颗粒粗糙度分析 % 输入: % xy - N×2 矩阵, 颗粒边界坐标, 第一列为x, 第二列为y % Nres - 等角度重采样点数, 建议 256 或 512 % K0 - 形状与粗糙度分界谐波阶次, 例如 8 或 12 % 输出: % stats - 结构体, 包含Ra, Rq, Rnorm, 振幅谱等 % % 依赖: 无额外工具箱 if nargin 2 || isempty(Nres), Nres 256; end if nargin 3 || isempty(K0), K0 8; end % 转到极坐标 cx mean(xy(:,1)); cy mean(xy(:,2)); dx xy(:,1) - cx; dy xy(:,2) - cy; r hypot(dx, dy); th atan2(dy, dx); % 统一角度到 [0, 2*pi) th(th 0) th(th 0) 2*pi; % 按角度排序 [th, idx] sort(th); r r(idx); % 等角度重采样, 闭合周期 thUniform linspace(0, 2*pi, Nres1).; thUniform(end) []; rInterp interp1(th, r, thUniform, linear); % 去均值 r0 mean(rInterp); rp rInterp - r0; % FFT 和单边幅值谱 Y fft(rp); N2 floor(Nres/2); Amp zeros(N21, 1); Amp(1) abs(Y(1)) / Nres; for k 2:N2 Amp(k) 2 * abs(Y(k)) / Nres; end % 奈奎斯特频率单独处理 if mod(Nres,2) 0 Amp(N21) abs(Y(N21)) / Nres; else Amp(N21) 2 * abs(Y(N21)) / Nres; end % 时域粗糙度指标 Ra mean(abs(rp)); Rq sqrt(mean(rp.^2)); % 频域粗糙度: K0以上谐波能量 kStart min(K0 1, N21); roughAmp sqrt(sum(Amp(kStart:end).^2)); Rnorm roughAmp / r0; % 打包结果 stats.r0 r0; stats.Ra Ra; stats.Rq Rq; stats.Rnorm Rnorm; stats.Amp Amp; stats.theta thUniform; stats.rInterp rInterp; stats.K0 K0; stats.Nres Nres; end这段代码的骨架看似简单但有几个细节值得展开。3.2 极坐标转换与重采样片段说明在极坐标转换部分我特意加了这一行代码th(th 0) th(th 0) 2*pi;目的很简单把角度统一到0到2π区间这样后续linspace(0, 2*pi, Nres1)与排序结果无缝对接。如果你不处理负角度排序后最前面的点可能是接近-π的而重采样网格从0开始中间就会留下一段拆断的间隙插值就会出问题。这是新手最容易踩的坑。插值部分用interp1的线性插值是因为颗粒轮廓点密度足够高时线性插值和样条插值结果差异很小而线性插值不容易产生过冲振荡。如果原始轮廓点很少或者颗粒表面有高频尖刺线性插值更稳定。另外重采样点数Nres的选择要匹配原始数据密度。我建议Nres取原始轮廓点数的四分之一到二分之一并限制在128到1024之间。512对我来说是最常用的默认值既能覆盖高谐波又不至于让高频段完全被噪声主导。3.3 FFT与单边幅值谱的细节实现FFT部分代码里那个循环其实就是把双边谱折叠成单边谱。有一个容易出错的边界情况当Nres为偶数时索引N21对应奈奎斯特频率分量的能量不需要加倍当Nres为奇数时没有真正的奈奎斯特频率最后一个有效分量仍需要加倍。我在代码里分别处理了这两种情况所以无论Nres取偶数还是奇数都能正确运行。实际测试中取偶数且为2的幂次是最省心的比如256、512、1024Matlab的FFT在这些长度上效率最高。另外振幅谱的索引含义要记住Amp(1)是直流分量Amp(2)是1阶谐波Amp(3)是2阶谐波以此类推。你也可以直接从Amp(3)开始看因为1阶谐波通常反映质心偏移物理意义有限。3.4 指标计算与输出结构最后一段计算Rnorm时我用K0以上所有谐波振幅的平方和再开方。K08的含义是第8阶以下的谐波算作形状能量第9阶及以上的谐波算作粗糙度能量。实际使用中我会同时输出完整振幅谱Amp方便后期做可视化。比如把Amp向量绘制成对数坐标的频谱图横轴为谐波阶次纵轴为振幅。这种图比一个单纯的Rq数字直观得多能一眼看出颗粒表面的主要波动落在哪个频段。输出结构体设计的重点是“不丢失中间信息”。即使当前只需要Ra也不要只返回Ra。因为你画出频谱后后续可能需要调整K0重新计算Rnorm如果没保存振幅谱就只能从头再算一次。代码的核心已经实现但还差一个验证用例。下面给出一段生成模拟颗粒并调用函数的脚本。% 模拟一个带宏观形状和微观粗糙度的颗粒轮廓 theta linspace(0, 2*pi, 2000).; rTrue 1.00 ... 0.10 * cos(2*theta - 0.5) ... 0.04 * sin(3*theta) ... 0.02 * cos(4*theta 0.8) ... 0.006 * cos(12*theta 1.0) ... 0.003 * sin(25*theta); x rTrue .* cos(theta); y rTrue .* sin(theta); % 人为抽取100个点, 模拟真实提取的不完整轮廓 idx round(linspace(1, 2000, 100)); xy [x(idx), y(idx)]; % 调用主函数 stats FourierRoughnessAnalysis(xy, 512, 8); % 显示关键结果 fprintf(等效半径 r0 %.4f\n, stats.r0); fprintf(Ra %.5f\n, stats.Ra); fprintf(Rq %.5f\n, stats.Rq); fprintf(归一化粗糙度 Rnorm %.4f\n, stats.Rnorm); % 绘制频谱 figure; semilogy(0:length(stats.Amp)-1, stats.Amp, o-); xlabel(谐波阶次 n); ylabel(振幅 (mm)); grid on;这段模拟数据里我把0.1的2阶分量、0.02的4阶分量当作“形状”把0.006的12阶分量和0.003的25阶分量当作“粗糙度”。如果K0取8算出来的Rnorm会比较好地呼应12阶和25阶分量的总能量如果K0误取4那形状里的4阶分量也会被算进粗糙度Rnorm会偏大。读者可以自己改一个K0跑一遍感受一下分界参数的影响。有个小技巧也顺便分享在模拟验证时尽量把各阶分量的数量级拉开。真实颗粒的形状分量比粗糙度分量常常高出30到100倍如果模拟数据里两者差异太小不利于观察计算方法是否正确。3.5 参数示例和验证用例验证代码跑完后理论上你会看到Ra和Rq主要由0.1和0.04这两项贡献Rnorm则由12阶和25阶这两项主导。我实际跑过的结果是Ra约为0.11Rq约为0.13Rnorm约等于0.007。如果算出来Rnorm在0.01以上多半是K0选小了或模拟数据里的“粗糙度”信号过强。判断实现是否正确还有一个更直接的频谱检查方法观察Amp(3)和Amp(5)这两个索引应该分别接近0.1和0.04因为2阶和4阶谐波对应这些位置Amp(13)和Amp(26)应该接近0.006和0.003。如果你的Amp符合这个规律说明FFT和单边谱换算没有出错。4. 常见坑从轮廓提取到频谱分析的故障排查4.1 角度排序和零度起点的断裂我先前提过的负角度问题很多人第一次跑通代码后并没有报错但画出来的频谱图高频段会有奇怪的尖峰群。原因通常是重采样网格跨越了角度断裂处插值插出了不存在的陡峭跳变。更常见的例子是轮廓本身从第三象限开始atan2输出的初始角度在-2.8左右排序后这一堆点被排到了序列最前面而重采样网格的起点是0这就意味着在网格起点附近会有一段区域几乎没有原始点覆盖interp1被迫在两个间距很远的点之间做插值生成了一段过大的跳跃。处理办法就是代码里那段角度统一加2π保证角度序列的首尾不会落在重采样区间两端之外。4.2 非均匀采样点直接做FFT的隐性错误这是新手最常犯的错误之一。不少人直接从bwboundaries拿到像素级边界点后不重采样直接对原始半径序列做FFT。感觉上好像也没错毕竟不管点间距均不均匀只要排列顺序是闭合的FFT就能给出某个结果。但这个结果物理意义很弱。原因在于FFT默认输入是等间隔采样信号。如果实际轮廓点角度间隔不均匀等间隔FFT相当于把一个非均匀采样信号强行当均匀信号处理高频谱会混入虚假分量。点间距差异越大虚假高频分量越强。有一些文献用“非均匀FFT”来处理这种点列但实际工程项目里重采样到均匀网格是最简单也最可控的方案。需要提醒的是重采样也会改变频谱细节所以报告的粗糙度指标一定要说明重采样点数Nres。同一个颗粒用Nres256和Nres1024算出来的高频段振幅有差异是正常现象不能认为代码有问题。4.3 图像噪声造成的高频伪尖峰从切片图像提取轮廓时边缘可能有像素级的毛刺和噪声这些毛刺对应的谐波阶次通常在50阶甚至100阶以上幅度不一定小。它们不仅会拉高Rnorm还会污染整个高频谱的形态。处理这个问题的顺序应该是在图像阶段先做轻微的高斯滤波或形态学闭运算把亚像素毛刺去除如果轮廓点已经提取出来了再用移动平均或小波方法做平滑但平滑尺度要控制得很小避免把真实粗糙度信息也抹掉。我在实际处理CT图像时常用3×3的中值滤波然后再提取边界效果比较稳定。另外要注意如果图像分辨率不足比如单个颗粒像素直径只有200个像素那么真正的“表面粗糙度”信息几乎不可分辨此时傅里叶展开得到的高频段更多反映的是像素离散化噪声。这种情况下强行计算高频粗糙度没有太大物理意义不如只提取前20阶谐波来分析形状特征。4.4 凹形颗粒的r(θ)多值问题这是极坐标方法的天花板。当颗粒形状明显内凹时从质心发出的某一条射线可能穿过颗粒边界两次这时r(θ)不是单值直接排序和插值会产生歧义。遇到这类情况我有两个建议。第一个建议是在前处理中判断颗粒凸度如果内凹深度超过等效半径的10%就考虑切换到复数形式的傅里叶描述子也就是把轮廓点的xiy当作一个复数序列来处理这样不再依赖r(θ)单值假设。第二个建议是针对具有轻微内凹但没有穿越质心的颗粒可以改用外接多边形重心替换质心通常能缓解一部分多值问题。值得一提的是真实砂土颗粒绝大部分是凸形或近似凸形。石英砂、玻璃珠、粉煤灰微珠这些材料极坐标半径方法几乎不会遇到凹形问题。真正棘手的是某些形状极不规则的珊瑚砂或人工破碎骨料这时要果断更换算法不要在极坐标上硬调。4.5 截止阶次与后续模拟对不上的问题在离散元模拟中如果接触模型需要输入一个综合粗糙度参数而K0选得不同标定出的摩擦系数可能差一倍。我遇到过一位同行用同一批颗粒轮廓数据一个课题组算出的粗糙度差两倍双方都认为自己的程序没问题后来一查发现是K0一个取了6一个取了20。这个问题没有标准答案但我建议这样做把K0当成模拟的一个输入变量不写成固定值。先做两组参考模拟一组K06一组K018看看宏观响应差异有多大。如果差异很小任何中间值都可以用如果差异显著就根据接触力学特征长度反推K0用颗粒尺寸和预期的微观起伏波长来定义分界。还有一种常用做法是定义截止波长而不是截止阶次。给定波长阈值λ0K0直接等于2πr0/λ0。这种表达方式更有物理直观也更容易在不同粒径的颗粒之间统一标准。5. 实操中的几个经验与建议最后分享几个我在实际项目中积累的体会。第一如果最近在写量化分析报告建议把振幅谱作为标配图放进去而不只是写Ra和Rq。频谱图能够在同一个坐标系里展现出颗粒的多尺度特征看图的人只要能理解横轴阶次对应波长就能获得比一个数字丰富得多的信息。不同批次样本对比时我通常直接在频谱图上叠加多条曲线粗糙度差异一目了然。第二重采样点数Nres和傅里叶截断阶次要一起记录放到论文或者报告的方法描述中。原始轮廓点数和重采样点数如果相差很多说明图像分辨率足够高计算结果可信度也高。如果原始边界点数本身还不到100个就别指望高谐波分量有意义那还不如直接用简单的椭圆拟合。第三关于程序性能Matlab的FFT对512点长度的序列处理几万个颗粒也只需要几秒到几十秒。瓶颈其实在于轮廓提取和插值。如果需要对一批颗粒批量计算我建议把坐标点数据预先保存成统一的二进制文件或mat文件避免每次重新读图。第四这个方法后续还可以扩展。比如对三维颗粒可以取多个正交剖面的轮廓分别做傅里叶分析再把各剖面的粗糙度指标做统计平均得到一个等效的三维粗糙度。虽然不能完全替代真正的表面形貌扫描但在常规颗粒材料研究中已经足够实用。计算岩土颗粒粗糙度这条路方法本身并不复杂复杂的是对边界条件和参数意义的理解。把傅里叶展开吃透配合稳定的轮廓提取流程比调了多少参数都有用。希望这份代码和踩坑记录能帮你少走几步弯路。