
简介这份资源面向光学测量、干涉检测与图像处理方向的学习者和研究人员聚焦相位提取与相位解包裹这一核心环节提供四步相移法程序与最小二乘法相位解包裹程序的配套实现。四步相移法通过多幅相移图求解相位主值最小二乘法解包裹则用于恢复连续相位分布两者组合可构成较完整的光学测量数据处理流程适合需要动手复现算法、验证方法效果的中高级读者参考。压缩包共7个文件以4个bmp图像文件、2个m脚本文件和1个db文件为主图像文件可作为算法输入样例m脚本承载相移与解包裹的核心计算逻辑整体约526KB体量轻便便于快速运行与调试。目前已有1337人学习下载说明该方案在同类资源中具有一定认可度。读者可据此理解相移法求解与最小二乘解包裹的衔接方式对照样例数据验证算法流程并在此基础上迁移到自己的干涉图或条纹图处理任务中。1. 四步相移法配最小二乘解包裹从条纹图到连续相位的完整链路拿到四张相位差固定为 π/2 的条纹图想恢复出连续、无跳变、能直接用于面形测量或干涉计量的相位分布这条链路的核心就是两件事四步相移法把相位主值算出来最小二乘法把被 arctan 折叠进 (-π, π] 的相位重新展开。标题里的“四步相移法程序”和“最小二乘法相位解包裹程序”不是两个孤立脚本而是一条流水线的上下游——前者负责逐像素求包裹相位后者负责在二维网格上把 2π 跳变抹平。做过干涉测量、数字全息、结构光三维重建的人对这套组合不会陌生它适合已经能采到相移图、但卡在“相位一展开就出现拉线或整片偏移”的从业者。下面按我实际写程序的顺序把公式、代码、参数和翻车点一次讲透。2. 四步相移法程序从四张条纹图算出包裹相位2.1 为什么选四步而不是三步或五步相移法家族里三步法只需要三张图采集快但对相移误差和背景光强波动很敏感五步法抗误差能力强可要多采两张图动态测量时反而容易引入运动伪影。四步法取相移量 0、π/2、π、3π/2在抗误差和采集成本之间取了一个工程上最常用的平衡点。它的强度表达式可以写成I_k(x,y) A(x,y) B(x,y)·cos[φ(x,y) (k-1)·π/2]k 1,2,3,4其中 A 是背景光强B 是条纹对比度φ 就是待求的包裹相位。把四式联立A 和 B 在相减过程中被消掉只剩 φ 的三角函数组合。这一步的选型理由很直接四步法对背景 A 和调制度 B 的慢变不敏感只要相移量准确逐像素解出的 φ 精度足够支撑后续解包裹。常见做法是先用四步法拿到主值再上最小二乘做全局展开而不是一上来就搞复杂的多波长或傅里叶变换。2.2 四步相移的逐像素公式与实现四步法的标准解相公式是φ_wrapped arctan2(I4 - I2, I1 - I3)这里用 arctan2 而不是普通 arctan是因为要保留象限信息否则解出的相位会丢一半。下面是我一般会写的 Python 实现输入是四张同尺寸的灰度图输出包裹相位矩阵。import numpy as np import cv2 def four_step_phase_shift(i1_path, i2_path, i3_path, i4_path): # 以灰度方式读入四张条纹图保证尺寸一致 i1 cv2.imread(i1_path, cv2.IMREAD_GRAYSCALE).astype(np.float64) i2 cv2.imread(i2_path, cv2.IMREAD_GRAYSCALE).astype(np.float64) i3 cv2.imread(i3_path, cv2.IMREAD_GRAYSCALE).astype(np.float64) i4 cv2.imread(i4_path, cv2.IMREAD_GRAYSCALE).astype(np.float64) # 分子分母分别对应 sin 和 cos 分量 numerator i4 - i2 # 对应 sin(phi) denominator i1 - i3 # 对应 cos(phi) # arctan2 输出范围 (-pi, pi]即包裹相位 wrapped np.arctan2(numerator, denominator) return wrapped逻辑说明分子 I4 - I2 提取的是 sin 分量分母 I1 - I3 提取的是 cos 分量两者一比再取 arctan2背景 A 被彻底消掉。参数上唯一要盯的是四张图的相移量必须严格是 0、π/2、π、3π/2如果实际相移器有偏差比如压电陶瓷非线性解出的相位会有周期性纹波。失败时先看 numerator 和 denominator 是否同时接近零——那说明该像素对比度太低属于无效点后面解包裹要把它掩掉。2.3 读图、归一化和无效点掩膜实际程序里我不会直接把原图丢进公式。先做两件事一是把四张图统一转成 float64避免 uint8 相减时下溢二是算一个调制度图 B sqrt(numerator² denominator²) / 2用它做无效点掩膜。调制度低于阈值的像素相位是噪声硬解出来只会污染整片解包裹结果。阈值一般取调制度最大值的 5% 到 10%具体看条纹对比度。这一步不做后面最小二乘会把噪声点当成真实跳变去拟合出现整片倾斜。modulation np.sqrt(numerator**2 denominator**2) / 2.0 mask modulation (0.05 * modulation.max()) # 无效点标记 wrapped[~mask] 0 # 无效点相位置零后续解包裹跳过参数说明0.05 这个系数不是固定的条纹质量差就调到 0.1宁可多掩掉一点也别让噪声点参与解包裹。掩膜之后建议把 wrapped 存成 .npy方便解包裹程序单独调试不用每次重读四张图。3. 最小二乘法相位解包裹程序把 2π 跳变抹成连续面3.1 最小二乘解包裹在解什么包裹相位 φ_wrapped 和真实相位 φ 的关系是 φ_wrapped φ - 2π·kk 是整数。解包裹的本质是给每个像素找一个整数 k让相邻像素的相位差尽量小。最小二乘法的思路不是逐行逐列去数跳变而是把“相邻像素展开后差值应等于包裹差值的主值”写成一个全局最小化问题最后归结为求解一个大型稀疏线性方程组。它的好处是对局部噪声和残差鲁棒不会像洪水填充那样一个点错就整片错代价是要解方程内存和计算量比路径跟踪法大。二维最小二乘解包裹的离散泊松方程形式是ρ(i,j) [Δx_w(i,j) - Δx_w(i-1,j)] [Δy_w(i,j) - Δy_w(i,j-1)]其中 Δx_w 和 Δy_w 是包裹相位在 x、y 方向的一阶差分再取主值。解这个泊松方程得到的就是最小二乘意义下的展开相位。3.2 用离散余弦变换快速解泊松方程直接解稀疏方程组在百万像素级图像上很慢工程上常用 DCT离散余弦变换把泊松方程对角化几步变换就能出结果。下面是我常用的实现输入是包裹相位和掩膜输出展开相位。import numpy as np from scipy.fftpack import dctn, idctn def ls_unwrap(wrapped, maskNone): if mask is None: mask np.ones_like(wrapped, dtypebool) # 计算 x、y 方向包裹差分的主值 dx np.zeros_like(wrapped) dy np.zeros_like(wrapped) dx[:, :-1] np.angle(np.exp(1j * (wrapped[:, 1:] - wrapped[:, :-1]))) dy[:-1, :] np.angle(np.exp(1j * (wrapped[1:, :] - wrapped[:-1, :]))) # 构造泊松方程右端项 rho rho np.zeros_like(wrapped) rho[:, 0] dx[:, 0] rho[:, 1:-1] dx[:, 1:-1] - dx[:, :-2] rho[:, -1] -dx[:, -2] rho[0, :] dy[0, :] rho[1:-1, :] dy[1:, :] - dy[:-1, :] rho[-1, :] -dy[-2, :] # 用 DCT 解泊松方程 rho_hat dctn(rho, type2, normortho) m, n wrapped.shape i np.arange(m)[:, None] j np.arange(n)[None, :] denom 2 * (np.cos(np.pi * i / m) np.cos(np.pi * j / n) - 2) denom[0, 0] 1 # 直流项单独处理避免除零 phi_hat rho_hat / denom phi_hat[0, 0] 0 phi idctn(phi_hat, type2, normortho) return phi逻辑说明dx、dy 用 exp(1j·Δ) 再取 angle等价于把差分折叠回 (-π, π]这是最小二乘解包裹的标准预处理。rho 是把 x、y 两个方向的差分散度加起来构成泊松方程右端。DCT 把方程变换到频域后每个频率分量独立相除denom 就是离散拉普拉斯算子在 DCT 域的特征值。参数上要注意 denom[0,0] 必须单独置 1 并把 phi_hat[0,0] 置零否则直流项除零会得到 NaN整幅图全废。掩膜目前只在预处理阶段用如果要严格处理无效区域需要把掩膜纳入加权泊松方程那是另一个量级的复杂度一般测量场景先掩掉无效点再插值即可。3.3 解包裹结果的验证与残差检查程序跑完不能直接信。我一般做两个检查一是把展开相位重新包裹回去和原始 wrapped 逐像素比残差应该接近零二是看展开相位沿 x、y 的梯度正常面形应该是平滑的如果出现密集的细条纹状残差说明有残差点没处理干净。rewrapped np.angle(np.exp(1j * phi)) residual np.angle(np.exp(1j * (rewrapped - wrapped))) print(最大重包裹残差:, np.abs(residual).max())如果最大残差明显大于零优先查三处相移量是否准确、调制度掩膜是否太宽松、差分主值那一步有没有把边界像素算错。残差检查是解包裹程序的黑匣子不看这一步后面拿去算面形就是玄学。4. 四步相移加最小二乘解包裹的避坑与排查4.1 相移量不准导致解包裹出现周期性斜坡现象展开相位整体平滑但叠加了一层和条纹周期一致的波纹面形看起来像被规律性压弯。原因压电陶瓷或相移器实际步进不是精确的 π/2四步法公式假设相移严格等间隔偏差会直接调制到相位主值里。解决先用已知平面或干涉仪标定相移器的实际步进量把偏差代入广义相移算法或者改用对相移误差不敏感的五步法。我一般会在程序里留一个相移误差补偿系数标定一次写死。4.2 调制度掩膜太松噪声点把解包裹带偏现象展开相位在低对比度区域出现整片倾斜或局部鼓包重包裹残差在这些区域明显偏大。原因调制度低的像素相位是纯噪声最小二乘是全局拟合几个坏点会通过泊松方程把误差扩散到邻域。解决把调制度阈值从 5% 提到 10% 甚至 15%先掩掉再解掩掉的区域最后用邻域插值补。血泪经验是宁可多掩别让噪声点参与全局方程。4.3 DCT 解泊松方程时直流项处理错误现象程序不报错但输出相位整体偏移一个常数或者直接出现 NaN。原因denom[0,0] 在离散拉普拉斯特征值里本来就是零直接相除会除零另外 phi_hat[0,0] 不置零反变换后整幅图会带一个无意义的常数偏置。解决denom[0,0] 置 1phi_hat[0,0] 置 0这两行必须成对出现。这个坑很隐蔽因为图像看起来“差不多对”只是绝对值不对做相对面形测量时容易蒙混过关做绝对测量就翻车。4.4 边界差分主值算错导致边缘拉线现象展开相位在图像最右列或最下行出现一条明显的拉线或突变。原因dx、dy 初始化成全零边界像素的差分没有正确填充泊松方程右端在边界处不闭合。解决dx 最后一列、dy 最后一行要按边界条件单独处理或者干脆把有效区域向内缩一个像素再解。我一般会在解包裹前把图像裁掉一圈边界代价是损失一行一列换来的是边缘干净。4.5 把包裹相位直接当展开相位用现象后续面形计算出现 2π 台阶或者相位解缠前后数值范围没变化。原因四步法输出的 arctan2 结果天然在 (-π, π]忘记调用解包裹程序或者解包裹程序因为掩膜全 False 直接返回了输入。解决在流水线里加一个断言检查展开相位的动态范围是否明显大于 2π如果还挤在 (-π, π] 里说明解包裹没生效。这个坑新手常踩因为包裹相位图看起来也有条纹容易误以为已经展开。5. 进阶技巧用质量图引导的加权最小二乘提升解包裹鲁棒性普通最小二乘对所有像素一视同仁但实际条纹图里总有低质量区域。一个我常用的进阶做法是引入质量图把调制度或相位导数方差作为权重解加权泊松方程。权重高的像素在方程里话语权大低质量区域被自然压制解包裹结果在噪声区不会乱跑。实现上不需要重写整个 DCT 框架可以把权重乘进差分项再用预处理共轭梯度法解加权方程如果嫌麻烦退一步的做法是先按质量图把低质量像素掩掉解完再插值效果也能提升一大截。验证加权效果的方法很直接构造一组带已知面形的仿真条纹图加不同信噪比的噪声分别跑普通最小二乘和加权最小二乘比对面形残差的 PV 和 RMS。我自己的习惯是每换一批测量对象先跑一遍仿真验证参数再上实测数据。实测里最值得盯的参数是调制度阈值和相移误差补偿系数这两个调好了四步相移加最小二乘解包裹这条链路在大多数干涉测量场景里都能稳定出连续相位。希望帮到你。本文还有配套的精品资源点击获取