简介光学干涉计量、频谱分析及成像技术中相位包裹常使连续相位被2π周期打断增大数据分析难度。针对这一难题资源包内以MATLAB实现为切入点通过建立相位梯度方程并采用最小二乘优化获得连续相位估计进而给出完整解包裹方案其中核心算法脚本通过构建连续相位模型并最小化包裹误差平方和来恢复相位连续性另一个示例脚本则用于生成模拟数据、调用算法并可视化验证可直观对比解包裹前后的相位曲线。整个压缩包共2个m文件均为MATLAB源码体积仅1KB可直接在MATLAB环境中打开运行代码精简便于分段调试目前已有314人学习下载。读者可借此掌握相位解包裹的建模思路同时通过对示例脚本的反复修改将这一方法迁移到干涉计量、光学检测等实际场景中对光学测量与信号处理方向的学生或工程师尤为实用。1. 相位包裹图上的偶数倍问题为什么非最小二乘不可把干涉仪或者数字全息重建出来的相位图放大看到的不是面型而是一圈一圈等高线样的条纹每一条线都代表相位折叠了 2π。反正切函数的主值被限定在 (−π, π] 区间真实面型一旦跨过这个边界就会产生锯齿状的包裹相位。于是相位展开的目标变得很直接给每个像素加一个合适的 2π 整数倍恢复出连续的真实相位。传统路径跟踪法沿着低噪声区域逐点积分遇到残差就得绕路噪声一复杂绕错一条边整片面型就出现一条“缝合线”局部解虽然是精确的全局却不可靠。这类问题的稳健解法是把整个展开过程当成一个最小二乘问题来解不追求某一条积分路径的自洽而是让所有像素误差的平方和全局最小。这篇内容面向正在处理干涉测量、数字全息和结构光投影数据的工程师会从数学推导、快速实现、加权处理一直讲到验证方法让你能独立写出可复现的展开程序。2. 最小二乘公式的法方程从包裹梯度到泊松方程2.1 包裹相位差不是直接相减而是折叠后相减设真实相位是 φ(x,y)观测到的包裹相位是 ψ(x,y)两者之间满足ψ W(φ) W(v) v − 2π·round(v/2π)当采样满足奈奎斯特条件相邻像素的真实相位差小于 π这时候连续面型对应的离散梯度信息可以完整保留在包裹差的折叠结果里。具体做法是用Δx(i,j) W( ψ(i1,j) − ψ(i,j) ) Δy(i,j) W( ψ(i,j1) − ψ(i,j) )来估计真实梯度。这里必须先折叠再使用否则原始差分在 π 边界处会得到接近 2π 的错误梯度后续积分会直接错一整圈。《今》小心的是现实数据不会这么理想噪声、阴影和相位不连续会令数值出现矛盾像素路径不同得到的结果也不同。最小二乘不复用单条积分路径而是构造目标函数E(φ) ∑ [ (φ(i1,j) − φ(i,j) − Δx(i,j))² (φ(i,j1) − φ(i,j) − Δy(i,j))² ]这里两个求和里的 Δx、Δy 已经由包裹差折叠得到关于 φ 的方程是纯二次的因此极小值点可以由解线性方程组唯一确定。2.2 最小二乘公式的法方程与离散拉普拉斯算子对 E 分别对每个像素求偏导并置零得到线性方程组(Dxᵀ Dx Dyᵀ Dy) φ Dxᵀ Δx Dyᵀ ΔyDx、Dy 是一阶差分矩阵Dxᵀ、Dyᵀ 是它们的伴随运算。左边的作用等价于五点离散拉普拉斯右边等价于把包裹梯度做离散散度。也就是说均匀权重的相位展开最小二乘解等价于求解带 Neumann 边界条件的泊松方程∇²φ ρ ∂φ/∂n 0边界处没有外部像素法向导数自然为零这是 Neumann 边界的由来。由于拉普拉斯算子存在常数零空间解的绝对水平不可确定实际处理中把解的均值归零或固定参考像素观测值。计算散度 ρ 是实现的第一步。下面这段代码不涉及快速算法只演示边界和梯度的处理骨架import numpy as np def wrap_pi(v): 把相位差折叠回 (-pi, pi] return np.arctan2(np.sin(v), np.cos(v)) def wrapped_gradient_divergence(psi): h, w psi.shape gx np.zeros((h, w)) gy np.zeros((h, w)) # 水平、垂直方向的包裹梯度 gx[:, :-1] wrap_pi(psi[:, 1:] - psi[:, :-1]) gy[:-1, :] wrap_pi(psi[1:, :] - psi[:-1, :]) # 对包裹梯度取散度边界上使用 Neumann 假设 rho np.zeros((h, w)) rho[:, 1:-1] gx[:, 1:-1] - gx[:, :-2] rho[:, 0] gx[:, 0] rho[:, -1] -gx[:, -2] rho[1:-1, :] gy[1:-1, :] - gy[:-2, :] rho[0, :] gy[0, :] rho[-1, :] -gy[-2, :] return rho这段代码里gx 是水平方向上的包裹梯度gy 是垂直方向上的包裹梯度。散度公式 rho gx(i,j) − gx(i,j−1) 加上 gy(i,j) − gy(i−1,j) 在边界处只保留内侧一项这正对应于 Neumann 边界条件。wrapped_gradient_divergence 得到的就是泊松方程右端 ρ后续所有解法都从这里出发。2.3 为什么加权版本的约束会不同均匀权重目标函数默认每个误差项同等重要。光学测量面型的噪声区域和低调制区域显然不满足这一条件坏像素的误差会通过最小二乘平均扩散到邻域形成一种低频“隆起”。因此真实数据往往需要给每条边分配权重目标函数变成加权平方和对应法方程左侧的拉普拉斯算子也将变成加权拉普拉斯。这一区别贯穿整个相位展开流程也是后面第 4 章独立成篇的原因。符号上的几个关键对象可以先统一起来符号含义尺寸ψ包裹相位观测H×W取值区间 (−π,π]Δx / Δy水平 / 垂直包裹梯度(H,W−1) / (H−1,W)ρ梯度散度泊松方程右端H×Wwx / wy水平 / 垂直边权重与 Δx / Δy 同尺寸3. 用 DCT 在整幅图像上求解均匀权重最小二乘展开3.1 为什么是 DCT 而不是直接求逆或 FFT均匀权重的最小二乘展开已经变成解一个大型稀疏泊松方程。直接构造稀疏矩阵再求逆在 512×512 图像上已经能感到明显的延迟更大尺寸时内存和时间的增长都不可接受。DCT 方案的关键在于当权重全部相同时离散拉普拉斯算子的特征向量恰好是余弦基因此可以在变换域直接完成对角化解。不使用 FFT 而用 DCT是因为 Neumann 边界条件要求图像边界外法向导数为零这正好对应离散余弦变换二期DCT-II的隐含边界扩展方式。在频率域中泊松方程的解退化成一次逐点除法。离散余弦变换本身存在快速算法整体复杂度为 O(H·W·log(H·W))对 2048×2048 大帧也足够快。3.2 最小可运行的 DCT 相位展开代码在得到第 2 章的散度 ρ 之后只需要再做一次正变换、一次逐点滤波、一次逆变换。下面是完整的均匀权重解import numpy as np from scipy.fft import dctn, idctn def unwrap_ls_dct(psi): 均匀权重最小二乘相位展开输入输出都是 HxW 数组 h, w psi.shape rho wrapped_gradient_divergence(psi) # 对右端项做二维 DCT-II v dctn(rho, type2, normortho) # 使用对应离散拉普拉斯频率分母 m np.arange(h)[:, None] n np.arange(w)[None, :] denom 2.0 * (np.cos(np.pi * n / w) np.cos(np.pi * m / h) - 2.0) # (0,0) 频率对应常数项分母为 0单独处理 out np.zeros_like(v) out[1:, :] v[1:, :] / denom[1:, :] out[0, 1:] v[0, 1:] / denom[0, 1:] # 强制全局常数归零不影响相对面型 out[0, 0] 0.0 return idctn(out, type2, normortho)分母 denom 来源于把五点拉普拉斯算子的特征值写成余弦形式。m0 且 n0 时 denom 为 0该频率分量对应全局偏置解相位展开问题时这个偏置没有物理意义因此直接置零。这里有一个容易踩的坑替换成 scipy.fftpack 老接口时dctn 的 norm 参数行为与 scipy.fft 不完全一样直接搬代码会得到整体缩放 2 倍或 4 倍的结果。只需要直观比较展开后相邻差和原始包裹梯度是否一致即可发现。uniform 权重下如果出现整体幅度偏差改一下规范化参数。3.3 边界尺寸和频率轴方向的参数取舍DCT 结果要求输入图像边缘是自然延拓这在实际干涉图中往往不成立因此解在图像四边会相对更柔和不像中心区域那样贴合梯度。常见做法是展开前对图像边缘做少量 padding或者在边缘区域用加权的迭代解法进行修正。这里对比一下不同解法的复杂度与适用场景解法时间复杂度内存特点适用场景DCT 均匀解O(H·W·log(HW))需要 2~3 个与图像同尺寸数组理想数据、完整图像、快速原型稀疏直接解中小尺寸可接受大尺寸急剧上升需稀疏矩阵存储加权、掩膜、小图验证迭代法 PCG/CG每次迭代 O(HW)次数与条件数有关只需算子运算不必显式存矩阵加权、带掩膜、大图3.4 一个必须留意的 DC 项陷阱在 3.2 代码中out[0, 0] 被置为零。如果你刚好把 denom[0,0] 设为一个很小的数再参与除法解会整体偏离一个巨大常数之后再用two pi判断时容易造成视觉误导。更好的做法是像示例代码那样从频域直接禁止该项进入逆变换。DCT 解法的最大限制是它只接受均匀权重。真实光学数据里如果存在坏像素、遮挡区DCT 仍然会把错误平滑到邻域。这就是下一章要解决的加权问题。4. 加权最小二乘给质量图加权的展开算法4.1 从目标函数到加权拉普拉斯方程给每条边分配权重以后目标函数变为E_w(φ) ∑ wx(i,j)·(φ(i1,j) − φ(i,j) − Δx(i,j))² ∑ wy(i,j)·(φ(i,j1) − φ(i,j) − Δy(i,j))²对 φ(i,j) 求导得到的法方程左侧是一个加权拉普拉斯它的结构不再是五点在常数权重上的均匀组合而是每条边的权重直接参与求和。相当于原点的 diag 由四个方向权重相加邻接项则是相应边权重的相反数。这不能再通过 DCT 一步求解需要把它送入共轭梯度或预先条件共轭梯度迭代。4.2 用 scipy.sparse 构建加权拉普拉斯并进行 CG 求解权重矩阵 wx、wy 应该与包裹梯度 Δx、Δy 的尺寸一致。构建稀疏矩阵时目标是让 A 的行和为 0且对角线自然包含边权累积。下面代码可以直接作为加权最小二乘求解器from scipy.sparse import coo_matrix, diags, identity from scipy.sparse.linalg import cg def solve_weighted_ls(psi, wx, wy, reg1e-8, maxiter2000): h, w psi.shape n h * w idx np.arange(n).reshape(h, w) # 水平边像素(i,j) 与 (i,j1) row_h idx[:, :-1].ravel() col_h idx[:, 1:].ravel() # 垂直边像素(i,j) 与 (i1,j) row_v idx[:-1, :].ravel() col_v idx[1:, :].ravel() r np.r_[row_h, row_v, col_h, col_v] c np.r_[col_h, col_v, row_h, row_v] val np.r_[-wx.ravel(), -wy.ravel(), -wx.ravel(), -wy.ravel()] A coo_matrix((val, (r, c)), shape(n, n)).tocsr() d A.sum(axis1) A A - diags(np.ravel(d)) reg * identity(n, formatcsr) # 右端权重乘包裹梯度并按边符号分布 gx wrap_pi(psi[:, 1:] - psi[:, :-1]) gy wrap_pi(psi[1:, :] - psi[:-1, :]) b np.zeros(n) np.add.at(b, row_h, wx.ravel() * gx.ravel()) np.add.at(b, col_h, -wx.ravel() * gx.ravel()) np.add.at(b, row_v, wy.ravel() * gy.ravel()) np.add.at(b, col_v, -wy.ravel() * gy.ravel()) phi, info cg(A, b, rtol1e-6, maxitermaxiter) if info ! 0: print(CG 未收敛检查权重或正则项) return phi.reshape(h, w)这里 A.sum(axis1) 给出的恰好是每条边上权重绝对值的总和用它构造对角线后矩阵行和为零这是拉普拉斯算子所必需的。reg 对应一个极小的正则项除了避免数值退化外也能让某些完全被掩膜隔开的小区域不至于进入奇异。np.add.at 用于按边累积右端项。如果使用普通同一个像素被多条边共享时后面的赋值会覆盖前面的贡献结果会遗失大量边信息。这是边索引方式最容易出错的地方。4.3 加速技巧以 DCT 解作为预条件子加权最小二乘的收敛速度受到权重分布影响质量图里如果存在大面积零权重矩阵条件数会明显变差。一种实用且可靠的方案是把均匀权重的 DCT 拉普拉斯逆算子当作预条件子在每次 CG 迭代内用unwrap_ls_dct配合仅对应均匀拉普拉斯的解来改善搜索方向。实现上可以把unwrap_ls_dct封装为LinearOperator再交由cg调用。实践中这比直接使用原始稀疏矩阵迭代快数倍尤其在 1024×1024 以上的图像中。不过要注意预条件子需要按当前残差输入而不是相位输入因此封装时应直接接收右端向量 rho内部执行 DCT 正变换、频域滤波、逆变换再返回结果。若把整个 unwrap_ls_dct 流程直接套用会因为多算一次包裹梯度散度而破坏迭代语义。4.4 质量图权重怎么取权重原则上应当反映该像素邻域相位信息的可靠程度。常见方案包括用干涉条纹调制度 m(x,y) 作为基础质量值再映射为边权重或者用相位梯度的倒数梯度越大权重越低。一个可用的实际映射是wx m(x,i)² · mask这里 mask 把阴影、遮挡等无效区域置 0。调制度在低噪声区域接近 1在噪声区域会明显下降平方操作能进一步拉开差距。如果图像本身是复数干涉图也可以直接用振幅或相干系数。需要避免的做法是把权重设成 0 和 1 的二值图再去解这会引入新的强不连续效果不如连续权重稳定。5. 残差点、正则项和边界:影响最小二乘展开质量的真实因素5.1 残差点为什么让全局解看起来“水波纹”残差点是包裹相位在 2×2 环路上不闭合的像素点环路包裹梯度和不再为 0而是 ±2π 的整数倍。对路径跟踪法来说残差点是绕路的障碍对最小二乘来说残差点意味着矩阵本身不再完全相容。最小二乘不会拒绝这些不一致而是把它们平均分摊到整个求解区域宏观表现就是大片的低频起伏就像水面波纹。5.2 先用残差图定位问题区域对展开前数据计算残差图可以快速看出哪些区域不能信任。一个向量化的残差检测实现如下def residue_map(psi): h, w psi.shape a wrap_pi(psi[:-1, :-1] - psi[:-1, 1:]) b wrap_pi(psi[:-1, 1:] - psi[1:, 1:]) c wrap_pi(psi[1:, 1:] - psi[1:, :-1]) d wrap_pi(psi[1:, :-1] - psi[:-1, :-1]) return np.round((a b c d) / (2 * np.pi))返回的非零位置就是残差点位置。若残差点密集且互相连接成带说明该区域噪声已经破坏了包裹梯度的可恢复性。这种情况下即使加权也仍然存在物理上的信息缺失需要在处理流程中提前标注。5.3 三个关键参数的作用和设置方向参数作用常见设置方向wx / wy 权重映射指数控制低质量区域对平方误差的惩罚权重调制度取平方到四次方噪声严重时取高次正则项 reg保证加权拉普拉斯矩阵可解避免全掩膜区域奇异取梯度典型量级的 1e−8 到 1e−6CG 迭代容差 rtol决定解收敛到什么程度1e−6 对 2π 量级相位足够迭代过多反而放大数值噪声权重映射指数是一个非常有用的调节旋钮。增大指数会让低质量区域边缘的约束变得更弱但同时也会让求解区域内部产生梯度不连续具体取值需要结合残差密度和结果中的跳变数量判断。正则项不宜过大否则会改变泊松方程本身使结果偏向于平滑常见误区是把它当成滤波强度来调。5.4 空洞和断开区域的权重策略图像里遇到圆形遮挡或阴影不能直接删掉这些像素否则矩阵的连通结构会被破坏。正确做法是保留像素位置但把与该像素相连的所有边权重清零同时通过 reg 给对应节点一个极小的对角钳制。这样的结果在空洞附近不会出现异常突变掩膜内的相位值会被插值到平滑边界上。如果这些区域本就不需要相位值也可以在求解后按 mask 把对应结果丢弃。从实际效果看加权最小二乘的用途不是替代路径跟踪而是给路径跟踪或区域增长算法提供一个全局一致的初值。很多工程代码里会先运行一次加权最小二乘再用残差图去做局部修正两种策略配合使用远比单一大面积平滑效果好。6. 验证技巧用模拟面型检验展开结果的合理边界6.1 构造已知面型做一个闭环自检生成一个连续光滑的模拟面型加上可控噪声包裹后调用前三章代码再与真实面对比是确认算法实现正确的最快路径。模拟面型可以选双高斯峰因为它在中心和边缘有不同梯度能覆盖包裹频率的主要形态。y, x np.mgrid[0:256, 0:256] truth 5.0 * np.exp(-((x - 128) ** 2 (y - 128) ** 2) / (2 * 40 ** 2)) noisy truth 0.03 * np.random.randn(256, 256) wrapped wrap_pi(noisy) unwrapped unwrap_ls_dct(wrapped) valid truth 0.1 rmse np.sqrt(np.mean((unwrapped[valid] - truth[valid]) ** 2)) jump np.abs(np.diff(unwrapped, axis0)[valid[:-1, :]]) print(rmse, np.sum(jump np.pi / 2))这里 valid 掩膜只统计真实信号高于噪声底座的区域避免边界上低幅值相位对误差统计的干扰。jump 统计的是展开结果里仍然存在的陡峭突变正常情况下几乎为 0。如果 jump 数很大说明从包裹到展开的某些环节仍然丢失了整圈信息。6.2 三个比 RMSE 更早暴露问题的指标RMSE 是一个宏观指标但它会被少数几个错误整圈平均掩盖。实际验证时建议优先看三个更敏感的指标展开结果中相邻像素差超过 π/2 的像素占比展开前后重建的包裹相位与原始包裹相位差残差图中的残差是否还在原位置。第三个指标尤其重要如果展开后的相位重新包裹后与输入不一致说明解存在明显的逻辑错误应先去检查梯度折叠算子而不是调试后面的求解器。6.3 把 DCT 解作为加权迭代初值加权最小二乘的 CG 迭代可以从零向量启动也可以从一个更好的初值启动。将 uniform 权重的 DCT 结果作为初值传给 cg 求解器能减少 30% 到 50% 的迭代次数同时不会改变收敛终点。实现时只需要把 unwrap_ls_dct 的结果展平后作为 x0 参数传入 cg其余保持不变。这个技巧对 4.2 节的代码是直接可用的对质量图权重比较平滑的数据尤其有效。本文还有配套的精品资源点击获取