DIC数字图像相关MATLAB主程序：从图像序列到应变应力场全流程

简介：这是一套基于MATLAB实现的数字图像相关（DIC）主程序资源，面向材料科学、工程力学与生物医学领域的学习者与研究人员，用于通过连续图像序列非接触式测量物体表面的位移、应变与应力。资源包共161个文件，包含82个m脚本、43张jpg与19张png图像、10个html说明页、5个fig图形文件及少量txt文档，压缩包约6.01MB，代码与示例图像、界面文件配套齐全。内容覆盖图像预处理、角点检测与特征配对、互相关匹配、位移场计算、应变应力求解及结果可视化等完整流程，脚本注释与文件结构便于逐步调试运行。目前已有4164人学习下载，适合初学者理解DIC算法原理，也方便进阶用户按实验条件修改优化代码，快速搭建自己的光学测量分析流程。

1. 从一叠试件照片到全场应变云图：DIC 主程序到底替你省掉了什么

手头有一组拉伸试验的序列照片，试件表面喷了黑白散斑，相机固定、光源稳定，几十到几百张图按时间顺序排好。你想知道的不是“试件断没断”，而是每一帧之间材料表面每个点往哪个方向动了多少、局部应变集中在哪里、应力是怎么随着位移场演化的。拿卡尺量、拿引伸计量，只能得到几个离散点的数据，而数字图像相关（DIC，Digital Image Correlation）要做的，是把整幅图像变成一个密集的位移场，再从这个位移场里算出应变场，最后结合材料本构关系推出应力分布。

这套 DIC MATLAB 主程序，解决的就是“从图像序列到位移、应变、应力”这条链路。它不负责采集图像，也不负责喷散斑，它负责的是最核心也最容易翻车的一步：把序列图像读进来，划好感兴趣区域（ROI），做子区匹配，算出整场的 U、V 位移，再通过差分或最小二乘拟合得到应变，最后按你给定的弹性模量和泊松比换算成应力。适合谁用？做材料力学实验的研究生、做结构变形监测的工程师、以及想用 MATLAB 快速验证 DIC 算法而不是从零写匹配代码的人。如果你手里已经有序列图，缺的是一个能跑通、能调参、能出云图的主程序，那这份资源就是冲着你来的。

2. 主程序怎么读图、划 ROI、做子区匹配：三个环节的参数逻辑

2.1 图像序列的读取与预处理：别让第一帧就埋雷

DIC 主程序的第一步不是算，而是把图读对。常见做法是用imread配合dir批量读取，但这里有个容易被忽略的点：文件名排序。如果相机导出的是img1.tif、img2.tif……img10.tif，直接按字符串排序会变成img1、img10、img2，序列顺序全乱，后面算出来的位移就是跳变的。我一般会在读取后强制按文件名中的数字重新排一次。

% 读取指定文件夹下的图像序列，并按文件名中的数字排序 imgFolder = 'D:\DIC_test\specimen_01'; % 图像所在文件夹 imgFiles = dir(fullfile(imgFolder, '*.tif')); % 按需改成 *.png / *.jpg % 提取文件名中的数字用于排序，避免 img10 排在 img2 前面 numIdx = zeros(length(imgFiles), 1); for k = 1:length(imgFiles) token = regexp(imgFiles(k).name, '\d+', 'match'); numIdx(k) = str2double(token{end}); end [~, sortOrder] = sort(numIdx); imgFiles = imgFiles(sortOrder); % 读取第一帧作为参考图，并转为灰度 refImg = imread(fullfile(imgFolder, imgFiles(1).name)); if size(refImg, 3) == 3 refImg = rgb2gray(refImg); end refImg = im2double(refImg);

这段代码的逻辑说明：dir拿到的是文件结构体数组，regexp提取文件名里最后一组数字作为排序依据，sort返回的sortOrder才是真正的时间顺序。参数上，imgFolder换成你自己的路径，文件扩展名按实际格式改。转灰度是因为 DIC 匹配通常在灰度图上做，彩色图反而增加计算量且不提升精度。im2double把像素值归一化到 0 到 1，后续做归一化互相关时数值更稳定。

预处理里还有一个常被跳过但影响很大的步骤：ROI 的划定。主程序一般会提供一个掩膜或者矩形区域，让你只对试件表面做计算，避免背景区域参与匹配产生伪相关。我一般会先用第一帧画一个矩形 ROI，把试件边缘裁掉几个像素，因为边缘子区容易超出图像边界导致匹配失败。

2.2 子区匹配的核心参数：子区大小、步长、搜索窗口

DIC 的位移计算本质是：在参考图上取一个子区（subset），在变形图上滑动一个搜索窗口，找归一化互相关（NCC）或平方差（SSD）最大的位置，那个位置偏移就是位移。主程序里通常暴露三个关键参数：子区大小subsetSize、计算步长stepSize、搜索半径searchRadius。

子区大小决定了空间分辨率和噪声抑制的平衡。子区太小，比如 11×11 像素，对散斑质量要求极高，噪声一大匹配就飘；子区太大，比如 61×61，位移场会被过度平滑，应变集中区域算不出来。常见做法是取 21 到 41 之间，散斑颗粒直径大约占 3 到 5 个像素时，子区取 31×31 是比较稳的起点。

步长决定输出位移场的密度。步长等于子区大小时，输出的是不重叠的网格点；步长小于子区大小，输出的是重叠的密集场。做应变分析时，我一般会把步长设成子区大小的 1/4 到 1/2，比如子区 31、步长 8，这样应变计算有足够的数据点做局部拟合。

搜索半径要覆盖你预估的最大位移。如果两帧之间试件移动了 20 个像素，搜索半径至少给 25 到 30。给太小，匹配直接跑到窗口边缘，位移被截断；给太大，计算量上去，而且可能匹配到错误位置。一个实用技巧是先用大搜索半径跑前两帧，看最大位移量级，再回头收紧。

% DIC 主程序核心参数设置 subsetSize = 31; % 子区边长，奇数，常用 21~41 stepSize = 8; % 计算步长，越小场越密，计算量越大 searchRadius = 30; % 搜索半径，需大于预估最大位移 nccThreshold = 0.85; % 相关系数阈值，低于此值视为无效点 % 在参考图上按步长生成计算网格点 [imgH, imgW] = size(refImg); halfSub = floor(subsetSize / 2); gridX = halfSub + 1 : stepSize : imgW - halfSub; gridY = halfSub + 1 : stepSize : imgH - halfSub; [gridX, gridY] = meshgrid(gridX, gridY); pointsRef = [gridX(:), gridY(:)]; % 每个点的参考坐标

参数说明：nccThreshold是过滤无效匹配的阀门，低于这个相关系数的点直接标记为 NaN，不参与后续应变计算。这个阈值不要设太低，0.7 以下基本是噪声在主导；也不要设太高，0.95 以上在散斑质量一般时会把大量有效点误杀。0.85 到 0.9 是常见区间。

2.3 从位移场到应变场：差分、拟合与应力换算

拿到每个网格点的 U、V 位移后，应变不是简单地对相邻点做差分就完事。位移场里含有匹配噪声，直接差分会把噪声放大成应变的剧烈波动。主程序里通常提供两种方式：一种是局部最小二乘拟合，在每一个点周围取一个小窗口，拟合一个平面，平面的斜率就是应变；另一种是全局平滑后再差分。

我一般会优先用局部最小二乘，窗口大小取 5×5 或 7×7 个计算点。窗口太小，噪声压不住；窗口太大，应变集中区域被抹平。这个窗口和子区大小是独立的，子区管匹配，拟合窗口管应变平滑。

应力换算相对直接：在弹性阶段，用胡克定律，平面应力状态下 σx = E/(1-ν²) (εx + ν εy)，σy 同理，τxy = G γxy。主程序里会让你输入弹性模量 E 和泊松比 ν，然后逐点算出应力分量。这里有个前提：这套换算只适用于线弹性、小变形。如果你的试件已经进入塑性或者大变形，应力结果只能当参考，不能当定量结论。

% 局部最小二乘拟合计算应变场（以 ex 为例） strainWindow = 5; % 拟合窗口半径，单位是计算点 E = 210e3; % 弹性模量，单位 MPa nu = 0.3; % 泊松比 h = stepSize; % 计算点间距，对应物理尺寸需按标定换算 % 假设 U 为位移场矩阵，尺寸与 pointsRef 网格一致 % 对每个有效点，取其邻域做平面拟合 [ny, nx] = size(U); ex = nan(ny, nx); ey = nan(ny, nx); gxy = nan(ny, nx); for i = 1 + strainWindow : ny - strainWindow for j = 1 + strainWindow : nx - strainWindow % 取邻域位移 uPatch = U(i-strainWindow:i+strainWindow, ... j-strainWindow:j+strainWindow); vPatch = V(i-strainWindow:i+strainWindow, ... j-strainWindow:j+strainWindow); if any(isnan(uPatch(:))) || any(isnan(vPatch(:))) continue; % 含无效点的邻域跳过 end % 平面拟合求偏导，此处用简单差分近似示意 [du_dx, ~] = gradient(uPatch, h); [~, dv_dy] = gradient(vPatch, h); [du_dy, ~] = gradient(uPatch, h); [dv_dx, ~] = gradient(vPatch, h); ex(i,j) = mean(du_dx(:)); ey(i,j) = mean(dv_dy(:)); gxy(i,j) = mean(du_dy(:) + dv_dx(:)); end end % 平面应力状态下的应力换算 sigmaX = E / (1 - nu^2) * (ex + nu * ey); sigmaY = E / (1 - nu^2) * (ey + nu * ex); tauXY = E / (2 * (1 + nu)) * gxy;

逻辑说明：gradient在这里做的是数值偏导，h是计算点之间的物理距离，需要你根据标定板或者试件已知尺寸换算。如果h用的是像素单位，那算出来的应变是“每像素”的应变，不是真实应变，这一步是很多人翻车的地方。应力换算的三个公式对应平面应力线弹性，E和nu按材料实际值填。算完之后，无效点仍然是 NaN，画云图时会自动留白。

3. 避坑与排查：DIC 主程序跑不出结果的五个血泪经验

3.1 位移场出现大面积 NaN，相关系数普遍偏低

现象：跑完主程序，位移云图上一大片空白，检查nccThreshold发现大量点低于阈值。原因通常不是算法错，而是散斑质量或者光照变化。散斑颗粒太大、太小、对比度不够，或者序列图像之间光源闪烁导致灰度整体漂移，都会让 NCC 峰值不明显。解决：先单独看第一帧和第十帧的直方图，如果均值差超过 10 个灰度级，说明光照在变，需要做灰度归一化或者重新采集。散斑方面，颗粒直径占 3 到 5 个像素是经验值，喷得太密或太稀都要重做。

3.2 位移场整体偏移一个固定值，像是“零点漂移”

现象：刚体平移试验里，试件没动，但算出来的 U 场有一个非零常数。原因多半是参考子区和变形子区的坐标对应关系搞错了，比如网格点生成时多加了半个子区，或者图像读取时第一帧不是真正的参考帧。解决：用已知刚体平移量做标定，比如手动把第二帧整体移 5 个像素，看算出来是不是 5。如果不是，检查pointsRef的生成逻辑和imread的帧顺序。

3.3 应变云图噪声极大，像雪花屏

现象：位移场看起来还算平滑，但应变场全是高频抖动。原因是位移场里的匹配噪声被差分放大了。解决：两条路，一是增大子区尺寸，从 31 加到 41 或 51，用空间换噪声；二是增大应变拟合窗口，从 5 加到 7 或 9。如果还不行，说明散斑质量本身不够，匹配噪声底子太大，再怎么平滑也是徒劳。

3.4 应力结果量级明显不对，差几个数量级

现象：算出来的应力是几百 MPa 或者零点几 MPa，和材料实际强度对不上。原因几乎都是单位问题。位移场如果是像素单位，应变就是无量纲的“像素应变”，但如果你在应力公式里用的E是 MPa，而h没有换算成毫米，那应力就全错了。解决：先确认标定比例，比如 1 像素对应多少毫米，把h换成物理单位，再检查E的单位是 MPa 还是 Pa。我一般会在主程序开头强制写一行单位注释，跑之前先读一遍。

3.5 计算速度慢到无法接受，一帧要跑几分钟

现象：网格点几千个，每个点都在全图搜索，MATLAB 循环扛不住。原因是搜索策略太暴力，没有用频域加速或者预测搜索。解决：常见做法是用normxcorr2做模板匹配，它底层是频域实现，比手写双重循环快很多；另一个技巧是用上一帧的位移作为下一帧的搜索中心预测，把搜索半径缩小到实际位移的 1.5 倍左右。如果还是慢，把步长从 8 调到 16，计算点数量直接降到四分之一。

4. 进阶用法：批量处理、位移场后处理与结果验证

4.1 批量跑完整个序列，自动保存每帧的位移和应变

单帧跑通之后，下一步是把整个序列串起来。主程序一般会提供一个循环框架，但这里有个细节：参考帧的选择。常见做法有两种，一种是固定第一帧为参考帧，后续每一帧都和第一帧匹配，适合小变形累积；另一种是相邻帧匹配，适合大变形但需要把增量位移累加。我一般会先用固定参考帧跑一遍，如果后面帧的相关系数掉得厉害，再切到增量模式。

% 批量处理整个序列，固定第一帧为参考 nFrames = length(imgFiles); U_all = cell(nFrames, 1); V_all = cell(nFrames, 1); for f = 1:nFrames curImg = imread(fullfile(imgFolder, imgFiles(f).name)); if size(curImg, 3) == 3 curImg = rgb2gray(curImg); end curImg = im2double(curImg); % 调用主程序的核心匹配函数，返回 U、V 场 [U, V, nccMap] = dic_match(refImg, curImg, ... pointsRef, subsetSize, searchRadius, nccThreshold); U_all{f} = U; V_all{f} = V; fprintf('Frame %d / %d done, valid points: %d\n', ... f, nFrames, sum(~isnan(U(:)))); end % 保存结果，避免重复计算 save('dic_results.mat', 'U_all', 'V_all', 'pointsRef', ... 'subsetSize', 'stepSize', 'searchRadius');

逻辑说明：dic_match是主程序里的核心函数，输入参考图、当前图、参考点坐标和三个参数，输出位移场和相关系数图。fprintf打印每帧的有效点数，方便你实时判断有没有帧跑崩。最后用save把结果存成 MAT 文件，后面调参重算时不用再跑一遍匹配。

4.2 用刚体平移和已知应变做验证，别直接信云图

DIC 结果看起来漂亮不代表对。我一般会做两个验证：第一，刚体平移验证，把一张图人为平移整数像素，跑 DIC，看位移场是不是均匀等于平移量，误差应该在 0.01 到 0.05 像素量级；第二，已知应变验证，用模拟散斑图施加一个已知的线性位移梯度，看算出来的应变和设定值差多少。这两个验证跑通，再去看真实试验的云图，心里才有底。

4.3 位移场后处理：平滑、插值和应力集中识别

原始位移场总带噪声，直接算应力会放大。我一般会在应变计算前做一次中值滤波或者高斯平滑，窗口取 3 到 5 个计算点。插值方面，如果网格点不够密，可以用griddata插到更细的网格再画云图，但注意插值不产生新信息，只是让图好看。应力集中识别靠的是应变梯度，我通常会把 von Mises 应力或者最大剪应变单独画出来，看集中区域是不是和试件几何缺口或者裂纹尖端对得上。对不上，先回头查 ROI 和标定，别急着下结论。

从那以后我每次跑 DIC 之前，都强制先用刚体平移图走一遍全流程，确认位移场没有系统性偏移、应变场没有异常噪声，才敢把真实试验的序列图丢进去。这个习惯帮我省掉了至少三次把错误结果写进报告的尴尬。希望帮到你。

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