四步相移解相位与MATLAB相位展开实战：从包裹相位到三维重建

简介面向光学测量、机器视觉与MATLAB算法研究人员这份资源提供了四步相移解相位的完整实现方案。内容基于四步相移干涉测量原理涵盖干涉图像获取、图像预处理、相位差计算与相位恢复等关键环节可帮助读者快速掌握最小二乘与傅里叶变换等解相方法适用于精密光学检测、微纳结构分析及生物医学成像等教学研究场景。压缩包内共10个文件以M脚本和BMP图像为主辅以TXT说明与ASV备份文件整体大小仅50KB轻量便于学习。已有536人浏览学习。资源包含可直接运行的MATLAB程序、四幅相移图像样本及简要说明文档读者可结合代码注释理解四步相移公式推导与实现细节并在此基础上扩展至三维全息显微等应用也便于进一步尝试迭代或机器学习类高级相位恢复算法。1. 四步相移解相位的起点四帧条纹为什么能算出相位做结构光三维测量或者干涉条纹分析的人看到“四步相移解相位”这个标题第一反应多半是相位主值用atan2就能算真正的麻烦在后面的“解”字。四步相移只需要四帧正弦条纹每帧之间相移 90 度用加减法和反正切把背景项消掉得到包裹在(-π, π]的相位主值但要接着把一圈一圈的相位跳变接成连续曲面才是工程里真正的分水岭。MATLAB 在这条路上的天然优势是矩阵运算直接对应图像从仿真验证到真实条纹处理十来行代码就能跑通。这篇文章按“数学原理 → MATLAB 实现 → 相位展开 → 多频外差 → 调试验证”的顺序把四步相移解相位这条链路完整讲清楚适合正在搭结构光原型、做干涉测量或者初次接触相位恢复的工程师。2. 四步相移的数学模型与相位主值求解2.1 四帧强度方程与差分组合四步相移的测量模型是投影仪或干涉仪产生一组正弦条纹相机采集到的第k帧强度为I_k(x, y) A(x, y) B(x, y) * cos(φ(x, y) δ_k)其中A是背景光强B是条纹调制度φ是待求相位δ_k是第k步相移量。四步相移取δ [0, π/2, π, 3π/2]展开后得到四组强度帧序相移量强度表达式I10A B cos φI2π/2A - B sin φI3πA - B cos φI43π/2A B sin φ把 I1 减 I3得到2B cos φ把 I4 减 I2得到2B sin φ。这两个差分项恰好构成一对正交分量相位主值就是φ atan2(I4 - I2, I1 - I3)这个式子的价值在于背景光强A被完整消掉乘性调制度B只在atan2里被短约不进相位。换句话说物体表面反照率不均匀、照明有渐变、条纹对比度不一样都不会直接影响相位结果。这是四步相移比三步相移更稳的核心原因。三步相移δ [0, 2π/3, 4π/3]少采集一帧但对相移器精度和探测器非线性更敏感四步多花一帧时间换来的是对背景和噪声更强的抑制。2.2 atan2 还是 atan相位主值的方向问题工程上最容易踩的坑是把atan2写成atan。atan的值域只有(-π/2, π/2)当真实相位落在第二、第三象限时atan((I4-I2)/(I1-I3))会把角度折叠回第一象限解出来的相位图出现对称性错误。atan2(Y, X)的值域是(-π, π]正好覆盖四步相移需要的全角度范围。在 MATLAB 中要特别注意参数顺序atan2(I4 - I2, I1 - I3)第一个参数是sin项第二个参数是cos项。如果写反相位会变成π/2 - φ的形式条纹方向看起来是对的但解出的高度场会整体反转。我的习惯是每次写完后用下面的仿真代码先跑一遍确认符号和值域都对再处理真实验数据。2.3 用 MATLAB 仿真确认符号和主值% 仿真参数图像大小、条纹频率、已知相位 H 256; W 256; [x, y] meshgrid(1:W, 1:H); f 8; % 沿 x 方向共 8 个条纹周期 phi_true 2*pi*f*x/W 0.3*sin(2*pi*y/H); A 100; B 70; % 背景光强和调制度 % 生成四步相移条纹 I zeros(H, W, 4); for k 0:3 I(:,:,k1) A B*cos(phi_true k*pi/2); end % 四步相移解相位主值 phi_rec atan2(I(:,:,4) - I(:,:,2), I(:,:,1) - I(:,:,3)); M hypot(I(:,:,4) - I(:,:,2), I(:,:,1) - I(:,:,3)) / 2; figure(Name, 四步相移仿真); subplot(1,2,1); imagesc(phi_rec); axis image; colorbar; title(包裹相位); subplot(1,2,2); imagesc(M); axis image; colorbar; title(调制度);代码里h 256; w 256;是为了看主值跳变足够清楚。f 8表示整个幅面横向 8 个条纹周期如果f太小相位展开没有难度验证不了 unwrap 的问题如果f太大单像素内条纹频率接近采样极限仿真结果会失真。phi_true同时包含 x 方向的线性相位和 y 方向的正弦微扰这样后面可以用它做误差对比。M是调制度图它表示条纹的局部对比度等于B的估计值数值越高说明该像素信噪比越好这个变量后面要用来当质量图。MATLAB 的atan2输出范围是(-π, π]所以phi_rec上能看到规则的锯齿状跳变这是包裹相位的正常形态。M图应该是渐变且与位置无关如果M出现明显的周期性条纹说明相移没有严格做到 90 度需要检查相移器或投影仪抖动。3. MATLAB图像处理实战从真实条纹图到相位主值3.1 读取四帧条纹并转双精度真实项目里相机输出的通常是uint8或uint16灰度图像。直接相减会在负数处截断成 0所以第一步必须是转double。我一般用im2double它会自动把 0-255 映射到 0-1相位计算对缩放不敏感所以这样处理不影响结果。files {frame1.tif, frame2.tif, frame3.tif, frame4.tif}; I zeros(size(imread(files{1}), 1), size(imread(files{1}), 2), 4); for n 1:4 I(:,:,n) im2double(imread(files{n})); end % 中值滤波去噪核选 3x3 I medfilt2(I, [3 3]);medfilt2属于图像处理工具箱。核大小推荐[3 3]再大会把条纹边缘的相位梯度抹平导致后续调制度偏低。如果没有图像处理工具箱可以用imgaussfilt(I, 1)代替但对椒盐噪声的抑制不如中值滤波。注意这一步要在四帧统一处理不要在每一帧上单独做不同的背景校正否则四帧的相对强度关系被破坏差分项也会被污染。常见的预处理参数如下表预处理项推荐参数说明中值滤波3x3去椒盐噪声保留条纹边缘高斯滤波sigma 1-2去随机噪声条纹密度高时选 1调制度掩膜取调制度直方图谷值去掉阴影、遮挡、离焦区域3.2 从真实条纹计算相位主值与调制度掩膜仿真里相位主值可以直接用atan2但真实场景有阴影、低反光区这些区域的B接近 0算出来的atan2值完全是噪声。所以在解相位前应当先计算调制度掩膜把无效像素设为NaN。I1 I(:,:,1); I2 I(:,:,2); I3 I(:,:,3); I4 I(:,:,4); num I4 - I2; % 2B*sin(phi) den I1 - I3; % 2B*cos(phi) phi atan2(num, den); M hypot(num, den) / 2; % 调制度阈值具体值要从直方图观察 mask M 0.08; phi(~mask) NaN;num和den是后面所有误差分析的关键。M值等于局部条纹调制度B不受背景光强影响。阈值0.08不是通用值它取决于投影仪亮度、相机曝光和物体反射率。我会先用histogram(M(:), 1000)看分布双峰之间的谷底就是阈值如果直接拍脑袋取 0.08可能在暗区留下大量低质量像素给相位展开埋雷。3.3 MATLAB画图检查相位跳变解完后不要急着做 unwrap先用imagesc和plot检查包裹相位。figure(Name, 四步相移检查); subplot(2,2,1); imshow(I1); title(I1 0°); subplot(2,2,2); imshow(I4); title(I4 270°); subplot(2,2,3); imagesc(phi); axis image; colorbar; title(包裹相位); subplot(2,2,4); imagesc(mask); axis image; colorbar; title(有效区域); figure; plot(1:W, phi(128,:)); ylim([-pi pi]); title(中间行包裹相位剖面);剖面线应当在-π和π之间形成锯齿状跳变。如果看到某段相位平滑超过π说明atan2的参数顺序写反或者相移方向跟预期不一致。如果剖面线上有大量毛刺说明调制度阈值太低噪声像素进入了解相位计算。4. 解相位不是一次atan2二维相位展开与质量图引导4.1 为什么不能在二维相位图上直接 unwrapMATLAB 提供的内置函数unwrap(p, [], dim)按指定维展开一维相位。二维图上最省事的写法是phi2d unwrap(phi, [], 2); % 先沿行展开 phi2d unwrap(phi2d, [], 1); % 再沿列展开这个写法在实验室理想条纹下可用但真实测量里会出问题。unwrap的判定逻辑是相邻点相位差超过π时加或减2π。如果有一小段噪声点把真实相位差推到了1.1πunwrap会误判成跳变并沿展开路径把2π误差一路传播下去。二维图先沿行再沿列等于把坏点误差传播到整个矩形区域最后相位图会出现从坏点向两边扩散的“条纹断层”。正确做法是引入质量图。调制度图M本身就是天然质量图调制度高的像素相位可信度高应该先展开调制度低的阴影和低反光区放在最后处理这样它们的误差不会污染大片区域。4.2 用调制度作质量图的简化解包裹实现我提供一个简化版的质量引导展开函数。它把所有像素按质量从高到低排序从最高质量的种子开始每处理一个像素都找它四邻域中质量最高的已展开像素用它估计当前像素的整数级次。function U unwrap_sorted(phi, Q) % 基于质量图排序的四邻域相位展开 % phi: 包裹相位, 无效像素设为 NaN % Q: 质量图调制度或自定义权重越大越可信 [H, W] size(phi); Q(~isfinite(phi)) -Inf; % 无效像素不参与展开 U phi; done false(H, W); [~, seed] max(Q(:)); done(seed) true; [~, order] sort(Q(:), descend); for t 1:numel(order) idx order(t); if done(idx) continue; end [x, y] ind2sub([H, W], idx); bestQ -Inf; bx 0; by 0; for ns [-1 0; 1 0; 0 -1; 0 1]. nx x ns(1); ny y ns(2); if nx 1 || nx H || ny 1 || ny W continue; end if ~done(nx, ny) continue; end if Q(nx, ny) bestQ bestQ Q(nx, ny); bx nx; by ny; end end if bx 0 continue; % 当前像素尚未连接到已展开区域留到下一轮 end % 用最高质量邻域修正 2pi 级次 U(x, y) U(bx, by) wrapToPi(phi(x, y) - phi(bx, by)); done(x, y) true; end end这个函数的核心是wrapToPi将当前像素与邻域像素的相位差压缩到(-π, π]然后叠加到已展开邻域值上。如果没有 Mapping Toolbox可以把wrapToPi(v)替换为mod(v pi, 2*pi) - pi。order按质量降序排列质量最高的种子先展开后续像素只在连接到已展开区域时才会被赋予级次。这个简化版用全局排序代替了真正的优先队列性能不是最优但正确性够用生产环境可以把排序改成优先队列思路完全一致。调用方式U unwrap_sorted(phi, M);需要留意函数里continue跳过的像素永远不会被再次尝试。如果图像有多个互不相连的独立区域孤立区域的高质量像素可能永远展开不了。实际测量中阴影断开的情况很多所以更稳妥的方式是对每个连通区域单独找到种子再对区域内部做质量引导展开。这里给的是最小可用版本理解了堆栈思路后可以自己扩展。4.3 展开结果的快速检查展开相位应当是连续的。可以用下面的差分检查gx diff(U, 1, 2); gy diff(U, 1, 1); figure; imagesc(abs(gx) abs(gy)); axis image; colorbar; title(展开相位梯度);梯度图上如果出现一条条亮线说明那里还残留2π级次误差。另一个更直观的方法是画剖面线plot(1:W, U(128,:));正常的展开相位剖面线应当平滑没有突发跳变。如果剖面线在某段出现陡峭的台阶优先去调制度图上查那个位置是不是落在了低质量区。5. 多频外差与三维重建四步相移的工程延伸5.1 从包裹相位到绝对相位的多频展开单频四步相移解出来的展开相位是相对的它差一个全局常数2πK。这个常数不随位置固定所以单频率无法直接得到绝对高度。工程上常见的做法是投影多组不同频率的条纹每个频率都做四步相移然后用低频相位展开高频相位。时间相位展开的分层递推公式是K_n round((f_n / f_{n-1} * Phi_{n-1} - phi_n) / (2*pi)) Phi_n phi_n 2*pi*K_n其中f_n是第 n 个频率的条纹周期数phi_n是它的包裹相位Phi_{n-1}是上一个已展开的低频绝对相位。MATLAB 实现如下freq [1, 8, 64]; % 由低到高三个频率 phi cell(1, 3); for n 1:3 % 假设已经采集并计算得到第 n 个频率的四步相移包裹相位 % 这里把 I1n~I4n 换成实际变量 phi{n} atan2(I4n{n} - I2n{n}, I1n{n} - I3n{n}); end Phi cell(1, 3); Phi{1} phi{1}; % 频率为 1全场只有一个周期不需要展开 for n 2:3 K round((freq(n) / freq(n-1) * Phi{n-1} - phi{n}) / (2*pi)); Phi{n} phi{n} 2*pi*K; end频率比freq(n)/freq(n-1)是关键参数。比值越大单次展开对相位噪声越敏感比值太小又需要更多频率层级。常见的组合是[1, 8, 64]或[1, 16, 32]。实际使用前我会在平面前拍摄一组静态条纹把每个环节的残留误差估算出来确保最相邻两级频率的相位噪声远小于π否则级次K会跳错。5.2 相位到高度映射与点云显示拿到绝对相位Phi{3}后三维重建最后一步是相位到高度映射。对常见的投影栅格系统相位差与高度近似成线性关系z a0 a1 * (Phi{3} / (2*pi)); % 显示三维散点 [X, Y] meshgrid(1:W, 1:H); scatter3(X(:), Y(:), z(:), 2, z(:), .);其中系数a0、a1需要对已知高度的平面做标定常见的做法是移动一个平面到多个高度对每个高度记录绝对相位然后用polyfit拟合线性系数。更精细的模型会加入镜头畸变和投影仪畸变但核心链路不变。5.3 四步相移在实际工程里的两个边界第一四帧采集有时间间隔如果物体在采集过程中移动四帧之间的相位关系被破坏解出来的相位会出现水波纹状误差。这种情况下只能用三步相移加更快的相机或者改用单帧条纹解调方法。第二调制度阈值不只是一个后处理开关它同时决定了测量范围。阈值太高会把低反光物体表面大面积挖掉阈值太低又让阴影噪声进入展开。我一般会用第 3 章的调制度直方图来定阈值并且在测试报告中同时记录阈值和有效像素占比。6. 验证四步相移解相位算法的三个MATLAB技巧6.1 用带噪声的仿真做端到端误差测试不要在相位域里加高斯噪声要在强度域里加这样更接近真实相机响应。构造已知相位生成四帧条纹后加入高斯噪声再过完整流程最终用wrapToPi对比误差。sigma 0.02; % 强度噪声标准差 In I sigma * randn(size(I)); phiN atan2(In(:,:,4) - In(:,:,2), In(:,:,1) - In(:,:,3)); MN hypot(In(:,:,4) - In(:,:,2), In(:,:,1) - In(:,:,3)) / 2; maskN MN 0.08; UN unwrap_sorted(phiN, MN); err wrapToPi(UN - phi_true); rms_err sqrt(mean(err(maskN).^2)); fprintf(RMS phase error: %.4f rad\n, rms_err);sigma要匹配实际 8bit 相机量化误差的量级一般取 0.01 到 0.03。如果 RMS 误差超过 0.1 rad就要先查maskN的阈值和四帧配准而不是怀疑算法。6.2 用调制度直方图定掩膜阈值不要在 mask 上拍脑袋。直接观察直方图双峰谷值histogram(M(:), 1000); xline(0.08, r-, Threshold);如果直方图只有一个峰说明条纹整体调制度偏低或者物体表面反光太强此时要先调整投影亮度或者相机曝光而不是压低阈值。6.3 把 NaN 处理放在 unwrap 之前调制度掩膜和NaN不是展开完再做必须在unwrap_sorted的输入阶段就处理。只要相位是NaN质量图里对应位置就要设成-Inf否则排序会把无效像素排在前面浪费大量时间而且可能从噪声区域错误引路。phi(~mask) NaN; U unwrap_sorted(phi, M); figure; imagesc(U); axis image; colorbar; title(最终展开相位);看到最终展开相位图上的条纹完全连接成连续曲面再回到调制度直方图确认一下阈值四步相移解相位这条链路才算真正跑通。本文还有配套的精品资源点击获取