简介面向实验力学、图像处理与应变测量领域研究者这份基于MATLAB的数字散斑生成文档提供了采用一阶与二阶形函数的两种完整算法。文档从位移函数定义出发详细说明散斑数量、大小、峰值强度及位移梯度等核心参数的设置方法并结合代码示例演示如何生成含非均匀线性位移的变形散斑图可应用于全息干涉、颗粒跟踪测速等实验场景。代码按输入参数、随机散布散斑、误差函数计算灰度、归一化输出图像等步骤组织便于理解从参数设定到图像生成的全过程二阶形函数部分还给出将数值位移改写为函数方程的思路以适配非线性应变的模拟需求。资源仅包含1个docx文件大小19KB内容紧凑、注释清晰适合正在开展数字图像相关方法研究或需要快速生成模拟散斑图像的本科生、研究生及工程技术人员。已有433人学习值得下载参考。1. 为什么仿真数字散斑里藏着变形梯度张量做数字图像相关DIC的人迟早会遇到同一个问题拿实验机拉出来的序列图验证匹配算法结果里永远混着加载偏心、光路畸变和传感器噪声误差说不清是算法带来的还是实验带来的。所以有经验的团队都会先做一步仿真——用MATLAB合成一对“真值已知”的变形散斑图位移场是写死在代码里的再拿自己的DIC程序反解回来误差一目了然。这份素材做的正是这件事用误差函数积分生成高斯型数字散斑分别按一阶、二阶形函数施加非均匀线性位移场。适合做DIC算法验证、应变测量参数标定以及刚接触散斑仿真的MATLAB图像处理初学者。散斑图看着只是黑白噪点但生成过程里对位移梯度的处理直接决定了后续应变计算的精度上限。2. 用误差函数合成初始散斑像素级强度如何由积分得到2.1 高斯散斑模型的像素强度计算散斑不是简单画几个白色圆点。单个散斑的光强分布通常用高斯函数近似圆心在随机位置 (X_s, Y_s)峰值强度 I0特征半径 a。像素 (i,j) 接收到的光强是这个像素面积内对高斯函数的积分。高斯函数没有初等原函数但积分结果可以用误差函数 erf 表示这就是素材代码里那一长串 erf 差分的来源。NX 512; NY 512; S 1200; a 4; I0 1; X NX * rand(S, 1); % S 个散斑中心的 x 坐标 Y NY * rand(S, 1); % S 个散斑中心的 y 坐标 [Xg, Yg] ndgrid(1:NX, 1:NY); % 像素坐标网格Xg 对应 iYg 对应 j I zeros(NX, NY); for s 1:S I I (erf((Xg - X(s)) ./ a) - erf((Xg 1 - X(s)) ./ a)) .* ... (erf((Yg - Y(s)) ./ a) - erf((Yg 1 - Y(s)) ./ a)); end I I0 * pi / 4 * a^2 * I; G mat2gray(I); imwrite(G, undeformed.tif);这段代码和素材里双重 for 循环的版本完全等价只是把“遍历像素”换成“遍历散斑”利用矩阵运算一次算完整幅图。erf((Xg-X(s))/a) 的意义是像素左边界到散斑中心的误差函数值减去右边界对应的值得到该像素覆盖散斑的面积比例。两个方向的面积比例相乘就是二维积分结果。外面的 pi/4a^2 是归一化系数保证单个散斑在全图的积分光强为 I0pi*a^2这样改变散斑数量 S 时图像平均亮度不会失控。2.2 输入参数与选型建议素材里通过 inputdlg 逐项弹窗输入图像尺寸、散斑数量、散斑大小和峰值强度。这个交互方式调试时很烦建议改成脚本开头的常量或函数参数。参数取值直接决定散斑质量常用的经验范围如下参数含义典型值对散斑图的影响NX, NY图像宽高像素512×512 或 1024×1024太小则子区分辨率不足太大则生成耗时剧增S散斑数量8002000过少则出现大片空白区域过多则散斑重叠严重a散斑特征半径像素35直接决定散斑平均直径影响DIC子区匹配精度I0散斑峰值强度0.51.0控制对比度过大会导致灰度饱和一个常见误区是把 a 设得很大想让散斑更“清楚”结果散斑连成一片自相关曲线变得平坦DIC 匹配时相关峰不尖锐。我一般先固定 a4、S1200生成后用第 6 章的自相关方法检查散斑平均尺寸再反过来微调参数。2.3 mat2gray 与 tif 输出的位深注意素材里用 mat2gray 把强度矩阵归一化到 01再 imwrite 存成 tif。这里有个暗坑imwrite 对 double 类型的 01 数据默认按 8bit 量化灰度只有 256 级。如果后续要做亚像素位移验证量化噪声会直接干扰相关峰定位。常用的做法是先转 uint16imwrite(uint16(G * 65535), undeformed_16bit.tif);注意不要加 imwrite 的压缩参数BMP 或未压缩 tif 在后续读取时速度更快。素材里直接写imwrite(G,undeformed.tif)在验证高精度位移场时建议改成 16bit 输出。matlab画图查看结果时用imagesc(G)加axis image就能直接看散斑分布是否均匀这一句在调试阶段比 imwrite 更常用。3. 位移场建模从数值位移到一阶形函数的算法原理3.1 数值位移与位移梯度的物理含义素材的主程序里有两个层面的一堆参数UX、UY 是整幅图像的刚体平移ux、uy、vx、vy 是位移梯度分量。刚体平移好理解梯度分量的意义要展开。设 u 方向位移场为 u(x,y)v 方向位移场为 v(x,y)一阶线性形式写作u(x,y) UX ux·x uy·yv(x,y) UY vx·x vy·y其中 ux 是 u 对 x 的偏导表示水平位移沿 x 方向的变化率uy 是 u 对 y 的偏导vx、vy 同理。这四个量正好构成位移梯度张量决定了变形图里每个像素的坐标映射。素材代码里生成变形图时用的坐标变换是x x - u(x,y) (1-ux)·x - uy·y - UXy y - v(x,y) -vx·x (1-vy)·y - UY写成矩阵形式就是素材里的 M [1-ux, -uy; -vx, 1-vy]。这个 M 矩阵描述了原始坐标到变形坐标的线性映射J det(M) 是映射的面积变化率。若 J0说明网格发生翻转生成的变形图会出现皱褶这在位移梯度较大时尤其要注意。3.2 形函数如何决定位移场的空间变化素材要求在 u 水平和 v 竖直方向的位移“以非均匀线性变化的方式进行位移”这里的关键是一阶形函数对应的位移场是位置坐标的线性函数位移值随位置变化所以位移场是“非均匀”的但位移梯度 ux、uy、vx、vy 是常量。换句话说整个图像区域的应变是均匀的。这在DIC子区匹配里是最常用的假设——子区内部的变形通常用一阶形函数描述。二阶形函数则更进一步位移场写成二次多项式u(x,y) a0 a1·x a2·y a3·x² a4·x·y a5·y²v(x,y) b0 b1·x b2·y b3·x² b4·x·y b5·y²此时位移梯度变成位置的函数应变的分布也不再均匀能表达弯曲、剪切带这类复杂变形。素材里说“对于其余 10 个参数只考虑位移为数值的形式”指的就是这些高阶系数在代码中作为常量输入不参与形函数方程的构造。3.3 从数值到位移函数方程的改造思路素材原始代码里位移是通过 inputdlg 输入的两个标量 UX、UY代入变形循环时直接参与坐标计算。要改成函数方程步骤是固定的先写出位移场的解析表达式再对 x、y 求偏导得到位移梯度场的解析式最后把循环里原来的标量运算替换成函数调用。例如要模拟一个水平方向线性变化的位移场UXf (x, y) 0.5 0.002 * x - 0.001 * y;对应的一阶梯度就是 ux0.002、uy-0.001。这样改的好处是位移场有了连续解析定义后续无论做亚像素插值还是和DIC反解结果对比都有明确的真值函数可用。4. 一阶形函数变形散斑实现与主循环改造4.1 把位移函数接进变形主循环基于第 2 章的初始散斑生成代码改造变形部分的循环。核心变化在 erf 差分的坐标左边界用当前坐标 (i,j) 的位移右边界用 (i1,j1) 的位移这一步对应素材里双重循环中 i1、j1 的处理保证像素覆盖范围在变形后仍然连续。% 一阶形函数位移场为位置线性函数非均匀线性变化 UXf (x, y) 0.50 0.002 * x - 0.001 * y; % u(x,y) UYf (x, y) -0.30 - 0.001 * x 0.003 * y; % v(x,y) % 解析偏导对应位移梯度常量 % ux 0.002, uy -0.001, vx -0.001, vy 0.003 Xc Xg - UXf(Xg, Yg); % 左/下边界变形坐标 Xc1 Xg 1 - UXf(Xg 1, Yg 1); % 右/上边界变形坐标 Yc Yg - UYf(Xg, Yg); Yc1 Yg 1 - UYf(Xg 1, Yg 1); Idef zeros(NX, NY); for s 1:S Idef Idef (erf((Xc - X(s)) ./ a) - erf((Xc1 - X(s)) ./ a)) .* ... (erf((Yc - Y(s)) ./ a) - erf((Yc1 - Y(s)) ./ a)); end Idef I0 * pi / 4 * a^2 * Idef; D mat2gray(Idef); imwrite(D, deformed.tif);这段代码里Xc 是变形后像素左边界对应的散斑图坐标Xc1 是右边界。位移场函数输入的是原始像素坐标输出的是该位置的水平位移两者相减得到变形后的坐标。对比素材原代码原来用i*uxj*uyUX硬编码的线性表达式被替换成了对 UXf 的函数调用因此 UXf 写成任何形式都不影响循环结构。注意这里位移函数的单位是像素负值表示向左、向上移动。4.2 强度守恒与 J 的幅值补偿素材代码里变形后没有对强度做幅值修正。严格来说变形造成局部面积变化散斑密度改变灰度强度应按面积比缩放。这个修正因子就是前面算的 Jdet(M)。当位移梯度只有 0.002 量级时J≈0.997影响很小可以忽略但位移梯度到 0.05 以上时强度偏差开始干扰相关计算。常见做法是把变形散斑强度乘上 1/abs(J)J (1 - ux) * (1 - vy) - (-uy) * (-vx); Idef Idef / abs(J);代码里 ux、uy、vx、vy 是位移场偏导在图像中心的取值。如果位移场是纯线性的一阶形函数J 全局为常量直接乘即可如果像第 5 章那样用二阶形函数J 随位置变化就需要逐像素计算 1/J 矩阵再点乘。我一般会保留这个修正虽然视觉上差别不大但DIC反解位移场时幅值不守恒会在相关曲面引入轻微的周期性误差。4.3 性能和边界条件的坑素材原始代码是三层循环嵌套遍历 i、遍历 j、sum 内部对 S 个散斑求和。512×512 图像配 1200 个散斑一次生成要跑几分钟。改成第 4.1 节的散斑循环后耗时降到秒级。如果散斑数量到 5000 以上可以用 parfor 替代 for s注意每个 worker 累加 Idef 时需要按列切片避免广播冲突。边界条件上有个常被忽略的点位移场会让散斑部分移出图像区域。位移为负时图像左侧会空出一条黑色带位移为正时右侧同理。素材代码没有处理这个问题直接 imwrite 的结果里会出现边缘黑边。常见做法是生成时把图像尺寸向外扩展若干像素生成完再裁剪回原始尺寸或者对超出边界的散斑直接跳过不参与累加。我在做 512×512、位移约 2 像素的仿真时会把计算网格扩成 516×516再截取中心 512×512黑边问题就消失了。5. 二阶形函数扩展多系数位移场的代码落法5.1 十二系数位移场与参数向量化二阶形函数完整的多项式归并后u、v 两个方向各 6 个系数共 12 个。素材代码只把初始散斑图和变形循环做了模板位移场部分需要自己补全。相比一阶形函数二阶改造的代码量没有增加多少核心还是把 UXf、UYf 的表达式换掉% 二阶形函数系数顺序为 a0,a1,a2,a3,a4,a5, b0,b1,b2,b3,b4,b5 p [0.2, 0.002, -0.001, 2e-6, -3e-6, 1e-6, ... -0.1, -0.001, 0.003, -1e-6, 4e-6, -2e-6]; UXf (x, y) p(1) p(2)*x p(3)*y p(4)*x.^2 p(5)*x.*y p(6)*y.^2; UYf (x, y) p(7) p(8)*x p(9)*y p(10)*x.^2 p(11)*x.*y p(12)*y.^2; Xc Xg - UXf(Xg, Yg); Xc1 Xg 1 - UXf(Xg 1, Yg 1); Yc Yg - UYf(Xg, Yg); Yc1 Yg 1 - UYf(Xg 1, Yg 1);和素材要求对应一阶形函数涉及的 ux、uy、vx、vy 以及两个平移量被归一化到 p 向量里参与函数方程计算其余高阶系数以数值形式放在 p 数组中素材中提到的“其余 10 个参数”走的就是这条路径不需要全部手动命名统一写成向量即可。这样处理的好处是想换一组系数时只改 p 数组不需要动下面的循环。5.2 系数矩阵的物理映射系数对应位移场物理含义典型量级a0, b0刚体平移整幅图像的常值位移03 像素a1, b1线性项u 方向沿 x 的拉伸/压缩1e-31e-2a2, b2线性项u 方向沿 y 的剪切效应1e-31e-2a3, b3二次项沿 x 方向的应变梯度1e-61e-5a4, b4交叉项剪切应变的空间变化1e-61e-5a5, b5二次项沿 y 方向的应变梯度1e-61e-5二次项系数取 1e-6 量级时在 512 像素的图像边缘产生的额外位移约 0.26 像素512²×1e-6既能让位移场出现可观测的非线性弯曲又不会让整像素偏移过大导致散斑匹配失败。如果要模拟更剧烈的局部变形把系数放大到 1e-5但此时要检查 J 是否出现负值一旦出现就必须减小系数。5.3 一阶与二阶选择过拟合与表征能力一阶形函数只有 4 个梯度参数反解时稳定但只能表达均匀应变二阶形函数能表达应变梯度适合弯曲试样、缺口件、剪切带这类变形高度局部化的场景。代价是反解参数从 4 个涨到 10 个以上相关搜索空间变大容易出现局部极值。做算法验证时我的习惯是先跑一阶确认流程通再换二阶如果 DIC 反解结果和设定位移场偏差大优先检查系数量级而不是怀疑算法。仿真阶段多花两分钟调系数比在实验数据上排查问题快得多。6. 生成后不自检等于白做散斑质量与DIC验证6.1 用自相关检查散斑平均直径生成完参考图先别急着做变形。用自相关函数检查散斑尺寸是否落在 35 像素区间这是DIC子区匹配最舒适的范围。F fft2(G - mean(G(:))); AC fftshift(ifft2(F .* conj(F))); AC AC / AC(NX/2 1, NY/2 1); profile_y AC(:, NY/2 1); half_width find(profile_y 0.5, 1, first) - 1; fprintf(散斑相关半径: %.1f px\n, half_width);AC 是归一化自相关图中心峰半高宽对应的像素数约等于散斑特征半径。这个值小于 2 说明散斑过细匹配时会受像素量化噪声影响大于 6 说明散斑过大相关峰太平亚像素插值精度下降。我一般会同时输出figure; plot(AC(NX/21, :))直接看剖面形状峰越尖锐越好。6.2 用互相关粗检整像素位移生成完变形图用整幅图像的互相关峰位置对比设定的平均位移能快速发现坐标变换方向是否有误。C fftshift(ifft2(fft2(G) .* conj(fft2(D)))); [~, idx] max(C(:)); [my, mx] ind2sub(size(C), idx); dx_est mx - NX/2 - 1; dy_est my - NY/2 - 1; disp([整像素位移估计: , num2str(dx_est), , , num2str(dy_est)]);设定的平均位移可以从位移场函数在图像中心处的值估算例如第 4 章的示例里中心位移约为 0.5 像素。如果互相关给出的位移方向相反或量级差很远多半是 Xc 和 Yc 的符号写反了。这个检查 30 秒就能跑完值得每次生成后都做。最后存图时统一用imwrite(uint16(D*65535), deformed_16bit.tif)保存变形图与参考图的 16bit 位深保持一致。读取后若发现边缘黑带优先检查图像边界是否预留了位移余量而不是改位移函数里的常数项。本文还有配套的精品资源点击获取