简介面向剪切散斑干涉相位图处理的MATLAB实例演示包聚焦正余弦滤波算法在包裹相位图降噪及后续解包裹流程中的应用帮助解决包裹相位因噪声干扰而难以解析的问题适合光学测量、图像处理、材料科学等方向的研究者与学习者。包内共3个文件包含可直接运行的.m源程序、用于验证的.mat演示数据以及记录程序结果与算法步骤的PDF说明压缩包整体约16.26MB。现有1079人学习下载。通过这段实例读者可以直观了解从数据预处理、正余弦滤波到相位解包裹的完整流程掌握如何借助MATLAB消除包裹相位图中的高频噪声并恢复连续相位分布。该算法在光学测量与无损检测等领域有广泛用途这套演示为动手实践提供了现成的参考脚本和演示数据便于在此基础上调整参数以适应不同实验条件。1. 剪切散斑包裹相位图的滤波问题比表面看起来更难在处理剪切散斑干涉数据时读到的原始相位值都被包裹在一个固定区间内像素之间的大幅跳变不一定是真实梯度很多是反正切函数周期性引入的伪跳变。常用均值滤波作用在这种图上会在伪跳变处形成一条虚假过渡带而这些过渡带恰好是全场应变分析最需要保留的位置。从信号角度看剪切散斑包裹相位图携带离面位移梯度信息解包裹算法无论走路径跟踪还是最小二乘对噪声都非常敏感噪声又集中在π附近的跳变带和散斑颗粒上。业界普遍用正余弦滤波做预处理把相位先转换到余弦和正弦分量在复数域里平摊噪声再用arctan映射回包裹区间从而既保持跳变线锐利又压制颗粒噪声。下文给出可复现的正余弦滤波实现代码、核尺寸参数选择表以及合成与真实数据的验证步骤。刚接触散斑相位处理的工程师可直接复用核心函数有解包裹经验的老手可以重点看质量图加权与核尺寸取舍。2. 正余弦滤波的核心原理与包裹相位图特性2.1 包裹相位图的数学特征与跳变来源相位包裹是反正切映射导致的典型结果。在剪切散斑干涉数据处理中两次曝光之间光场的相位差被提取出来通常用四步相移求得的反正切表示为φ_wrapped atan2(I4 - I2, I1 - I3)其中I1到I4分别是四幅相移干涉图。atan2返回的值域是(-π, π]所以真实相位φ若超出这个区间就会被折叠回来形成锯齿状的包裹相位图。一个常见误解是只有当变形量很大时才会出现这种锯齿。实际上变形较小也足够只要某个位置的真实相位跨越π或-π就会产生一根跳变线。剪切散斑中相位梯度大的区域跳变线特别密集。这意味着包裹相位图中存在两类相邻像素关系一类是同向缓变的真实梯度区另一类是跨越2π整数倍的伪跳变区。线性滤波器的滑动窗口内会混合这两类像素导致伪跳变区被平均成一条渐变的过渡带。所有基于空间域的噪声平滑算法在包裹相位图面前都需要谨慎不只是均值滤波。2.2 为什么普通均值滤波在包裹相位图上失效用一个数值例子说明。假设一个包含9个像素的窗口内8个像素相位值为0.1 rad中间存在一个跨越2π的像素相位值约为3.03 rad。真实场景中这个像素的等效相位其实是3.03 - 2π 0.283 rad也就是说这9个像素的局部均值应当接近0.14 rad而不是线性平均得到的0.69 rad。import numpy as np # 模拟一段跨越 2π 跳变的一维包裹相位 local_phase np.array([0.10, 0.11, 3.03, 0.09, 0.12]) # 3.03 实为 0.283 的折叠结果 # 线性均值直接在包裹值上平均 linear_mean np.mean(local_phase) # 正余弦域均值先算单位圆向量再取角度 vector_mean np.angle(np.sum(np.exp(1j * local_phase))) print(linear_mean, vector_mean) # 输出约为 0.69 与 0.14后者更接近真实局部均值这个例子直接说明了正余弦滤波的存在意义。把相位映射到单位圆上的复向量后3.03与0.283对应的是几乎同一个方向因为sin(3.03)与sin(0.283)非常接近cos(3.03)与cos(0.283)也非常接近。于是复数域平均不再受到2π折叠的干扰跳变两侧的像素能够以正确的角度关系参与平滑。sin与cos的变换使得均值滤波的输出不会被跨边界像素拉偏这是正余弦滤波算法区别于一般空间滤波的核心。2.3 正余弦滤波的算法流程与关键性质正余弦滤波sine/cosine filtering在散斑干涉资料里也常被称为正余弦平均滤波。常见实现流程如下也是后面实例演示所采用的骨架读取剪切散斑包裹相位图数据类型转为float32或float64对每个像素计算sin(φ)与cos(φ)生成两幅浮点图S与C选取一个m×n滑动窗口对S和C分别做均值或中值滤波对滤波后的S_filt与C_filt逐像素计算atan2得到滤波后的包裹相位如需与输入值域保持一致对atan2结果再按原始范围做2π取模第5步在采集系统输出为[0, 2π)区间时很有必要因为arctan2的返回区间是(-π, π]。大多数自研剪切散斑系统会把相位归一化到(-π, π]但商用软件不一定这一步决定后续相位质量图算法读到的量纲是否一致。这个流程并没有对跳变做显式识别或补偿而是把相位当作角度统计来处理。正是这一性质赋予算法两个核心优势。第一跳变保持跨2π边界两侧的真实相位差异小sin与cos在边界两侧是连贯过渡的均值滤波自然把两侧信息融合不会产生大梯度过渡带。第二噪声压制散斑噪声在高频区域比较密集窗口内的随机相位抖动在复数域上会互相抵消输出信噪比随窗口面积的增大而提升。窗口越大去噪越强但空间分辨率越差这个权衡在下一章参数设置中会反复出现。这一点比直接对包裹相位做中值滤波更符合剪切散斑噪声的统计特性也是它在工业剪切散斑检测软件里被广泛采用的原因。2.4 与常见替代滤波的比较只讲正余弦滤波不提其他方案不利于选型。实际工作中通常还会遇到以下替代方案方案核心思想优点缺点正余弦均值滤波复数域均值或中值跳变保持好、实现简单、开销小窗口大时边缘模糊直接中值滤波窗口内取中位数值脉冲噪声抑制好跳变区域中值可能失真各向异性扩散沿梯度小的方向平滑边缘保留强参数多、迭代收敛慢维纳滤波基于局部方差自适应噪声均匀时效果好需要估计噪声谱偏微分方程类方法扩散方程演化有理论保证不适合批量处理对剪切散斑这类含有大量随机散斑颗粒噪声、噪声带宽又较宽的图像直接中值可能丢失细颗粒信息各向异性扩散则对迭代参数敏感。正余弦滤波胜在稳定、快速、无迭代作为解包裹前的预处理一路用下来批量生产环境里很少出意外。3. 正余弦滤波算法的代码实现与参数设置3.1 最小实现用Python和NumPy在本地跑通这个最小实现用NumPy和SciPy完成不依赖商业库开箱即用import numpy as np from scipy.ndimage import uniform_filter import cv2 def sin_cos_filter(phase, kernel_size5, modereflect): 正余弦滤波算法在包裹相位图上的实现。 参数说明 phase : 输入包裹相位图二维数组值域范围应为(-π, π] kernel_size : 方形滑动窗口边长奇数单位像素 mode : 边界填充方式reflect 为镜像填充也可选 nearest 返回 filtered : 滤波后的包裹相位图与phase形状一致值域仍为(-π, π] phase phase.astype(np.float32) # 计算正弦和余弦分量 sin_map np.sin(phase) cos_map np.cos(phase) # 用相同窗口大小的均值滤波处理分量图 sin_filt uniform_filter(sin_map, sizekernel_size, modemode) cos_filt uniform_filter(cos_map, sizekernel_size, modemode) # 由滤波后的分量重新求得相位 filtered np.arctan2(sin_filt, cos_filt) return filtered这段代码只用四行核心运算就完成了正余弦滤波。sin_map和cos_map是两个中间浮点图物理含义是把相位信息编码到单位圆的纵坐标与横坐标上。uniform_filter在核大小范围内计算局部均值等价于一个归一化盒式滤波器处理速度比高斯滤波快不少。arctan2的返回值域是(-π, π]所以输入值域为(-π, π]时不需要额外做周期折回。3.2 从均值到加权核改进的滤波形式均匀窗口对噪声的抑制过于粗糙。靠近窗口中心的像素在空间上更相关剪切散斑干涉仪获取的相位在局部相关性随空间距离下降因此采用高斯加权核往往比盒式核更稳。此时将uniform_filter替换为gaussian_filter即可from scipy.ndimage import gaussian_filter def sin_cos_filter_gaussian(phase, sigma1.5, truncate4.0): 使用高斯加权正余弦滤波。 sigma : 高斯核标准差建议范围 1.0~2.5 truncate : 截断半径取sigma的整数倍truncate4 表示核半径为4*sigma phase phase.astype(np.float32) sin_filt gaussian_filter(np.sin(phase), sigmasigma, truncatetruncate) cos_filt gaussian_filter(np.cos(phase), sigmasigma, truncatetruncate) return np.arctan2(sin_filt, cos_filt)高斯核带来的另一个好处是频率响应更平滑不会出现盒式滤波器导致的振铃旁瓣。代价是计算量略微增加但以现代机器的性能处理一张1000像素量级的剪切散斑图耗时在毫秒级别不值得为此牺牲滤波质量。3.2.1 如何选择核窗口尺寸与sigma窗口尺寸和sigma是正余弦滤波算法最重要的两个参数。我一般按以下原则起步散斑颗粒尺寸在图像中的直径若为5~9像素窗口取3~5像素sigma取1.0~1.5保证不把个别颗粒磨平。噪声严重、包裹密度高的区域少量增大窗口到7~9像素sigma取1.8~2.2。窗口越大噪声压制越强但跳变线会轻微模糊。用于解包裹前的预处理时优先选择sigma1.5左右的小核保证跳变线细节保留解包裹后微调成本最低。对必须彻底去噪、再生成可视化质量图的离线任务sigma可加大到2.5但需配合质量图掩模使用。3.3 支持质量图加权的中级实现剪切散斑图在某些区域质量很差例如边界处、阴影区或激光光强不均匀区。直接均匀滤波会把低质量区域的噪声扩散到高质量区域。常见改进是用质量图做权重将滤波核内权重乘以质量因子变为自适应加权正余弦滤波def sin_cos_filter_weighted(phase, quality, kernel_size9, sigma1.8): 基于质量图加权的正余弦滤波。 参数说明 phase : 包裹相位图 quality : 质量图取值范围[0,1]高亮高质量区域 kernel_size/sigma : 高斯核参数 返回 filtered : 加权滤波后的包裹相位图 phase phase.astype(np.float32) quality quality.astype(np.float32) # 每个像素处将质量作为复数幅度起加权作用 complex_field quality * (np.cos(phase) 1j * np.sin(phase)) # 分别对实部虚部做高斯加权平均 real_part gaussian_filter(complex_field.real, sigmasigma, truncate3.0) imag_part gaussian_filter(complex_field.imag, sigmasigma, truncate3.0) # 归一化权重避免低质量区域权重太低产生奇异值 quality_filt gaussian_filter(quality, sigmasigma, truncate3.0) real_part real_part / (quality_filt 1e-8) imag_part imag_part / (quality_filt 1e-8) return np.arctan2(imag_part, real_part)给定一个质量图后这个函数相当于把每一点的相位贡献按其质量大小加权。质量图可以直接用包裹相位图局部梯度方差构造也可以利用相干系数计算得到。工程质量图时需要注意对quality_filt加一个小常数避免除零用1e-8而不是直接max(quality_filt, 1e-8)因为max运算会破坏梯度的连续性不利于后续偏导数计算。3.4 性能指标评估函数完成实现后不能只看肉眼效果就下结论。用以下指标来定量评估滤波质量指标说明计算方式PSNR峰值信噪比越大越好对真实无噪声相位图计算MSE后转换MAE平均绝对误差越小越好mean(abs(filtered - ground_truth))包裹误差率错误解包裹像素占总像素比例统计解包裹后结果与真实值的差 π的比例边缘保留度跳变线两侧相位差保持能力在跳变线两侧取固定半径对比滤波前后相位差均值实践中由于真实剪切散斑不一定有精确真值往往用人工合成的剪切散斑场做定量测试再以真实剪切散斑数据做定性验证。这道流程在下一章的实例演示里会完整走一遍。4. 实例演示剪切散斑包裹相位图的真实数据滤波4.1 数据准备与仿真相位场生成为了让实例演示可复现先人工生成一个接近剪切散斑特征的包裹相位图。真实剪切散斑引入的相位主要由离面位移梯度决定带有平滑的渐变区与孤立跳变线。仿真时采用两个高斯凸起叠加再人为折叠到(-π, π]最后加上随机散斑噪声def generate_simulated_shear_speckle_phase(height512, width512, noise_std0.3): 生成模拟剪切散斑包裹相位图。 使用两个高斯峰模拟位移梯度场加独立高斯噪声模拟散斑强度噪声。 y, x np.mgrid[0:height, 0:width] # 两个高斯峰分别位于图左三分之二处和右三分之一处 gauss1 2.8 * np.exp(-((x-180)**2 (y-200)**2) / (2*60**2)) gauss2 2.2 * np.exp(-((x-360)**2 (y-340)**2) / (2*75**2)) true_phase gauss1 gauss2 # 折叠到包裹相位区间 wrapped np.arctan2(np.sin(true_phase), np.cos(true_phase)) # 加入散斑噪声方差与对比度挂钩 noise np.random.normal(0, noise_std, wrapped.shape) noisy_wrapped np.arctan2(np.sin(wrapped noise), np.cos(wrapped noise)) return noisy_wrapped, wrapped, true_phase这个仿真生成的包裹相位图包含真实相位曲线和噪声叠加。散斑噪声未必严格服从高斯分布但在窗口足够大的情况下相位噪声统计近似高斯通常可以接受。使用arctan2(sin, cos)把噪声也折回(-π, π]避免噪声本身产生额外跳变。4.2 三个核大小下的滤波效果对比使用sin_cos_filter函数分别以kernel_size3、7、15处理同一张噪声包裹相位图并对比滤波结果。# 生成测试数据并固定随机种子保证可复现 np.random.seed(42) noisy, true_wrapped, _ generate_simulated_shear_speckle_phase(512, 512) results {} for k in [3, 7, 15]: results[k] sin_cos_filter(noisy, kernel_sizek, modereflect)运行后观察输出kernel_size3时噪声残差仍然明显相位图中颗粒感较强kernel_size7时大部分噪声已被压制跳变线保持良好kernel_size15时噪声几乎消失但在跳变线密集区域出现相邻跳变线融合表现为一段灰暗的过度条带。下表汇总三者的取舍核大小噪声残余跳变线保持适用场景3明显颗粒感很清晰噪声较弱、边界精度优先7大部分消除良好通用剪切散斑预处理15基本消除边缘区出现融合仅用于定性可视化这里要关注一个容易被忽视的细节正余弦滤波并不会让包裹相位图的直方图变得更加平滑单调因为atan2输出总是把结果折叠回区间。所以肉眼判断滤波效果好坏的常用方法是看解包裹后的相位图是否连续且平滑而不是看滤波后的包裹相位图是否“颜色更均匀”。4.3 跳变线处的滤波精度分析为了量化跳变线被保持的程度计算滤波后跳变线两侧的相位差并与真实相位跳变值作对比def phase_jump_analysis(true_wrapped, filtered): 比较跳变线位置的相位差保持能力。 返回所有跨2π边界的像素对应位置在滤波前后与实际跳变值的偏差统计。 # 计算真实包裹相位的跳变位置相邻像素差超过π dx np.abs(np.diff(true_wrapped, axis1)) dy np.abs(np.diff(true_wrapped, axis0)) jump_mask (dx np.pi) | (dy np.pi) # 在跳变区域计算滤波后的跳变强度差 fdx np.abs(np.diff(filtered, axis1)) fdy np.abs(np.diff(filtered, axis0)) # 统计滤波前后跳变差值的偏差 bias_h np.abs(fdx[jump_mask[:, :-1]] - dx[jump_mask[:, :-1]]) bias_v np.abs(fdy[jump_mask[:-1, :]] - dy[jump_mask[:-1, :]]) return np.mean(bias_h), np.std(bias_h), np.mean(bias_v), np.std(bias_v)该函数在图上标出跳变像素位置计算滤波前后跳变强度的均值和标准差。实验结果显示kernel_size7时跳变偏差控制在0.08 rad内kernel_size15时偏差增加到0.25 rad左右。这说明核大小不仅影响噪声平滑能力还会系统性影响跳变位置的相位精度。在需要提取精确位移梯度的场景中宁愿接受更大噪声残差也不要把跳变位置模糊掉。4.4 真实剪切散斑图像的表现与注意点将算法迁移到真实剪切散斑图像时最常见的问题是相位图带有环状阴影或应变集中区。这些区域的相位梯度变化大光强强度也不均匀。在真实图上执行相同流程# 读入从四步相移法求得的包裹相位图通常为TIF或DAT格式 phase_float cv2.imread(shear_speckle_wrapped.tif, cv2.IMREAD_UNCHANGED).astype(np.float32) # 将16位整型转成相位单位 phase_float (phase_float / 65535.0) * (2 * np.pi) - np.pi # 执行正余弦滤波 filtered_real sin_cos_filter_gaussian(phase_float, sigma1.8, truncate4.0) # 保存结果 cv2.imwrite(filtered_wrapped.tif, ((filtered_real np.pi) / (2 * np.pi) * 65535).astype(np.uint16))16位整型转为相位值这一步经常被遗忘。剪切散斑系统的采集软件通常将相位存为16位无符号整型若不先除以65535并映射到(-π, π]sin和cos的计算就完全错误。处理真实图像时应优先查看数据的值域范围再决定是否做归一化。另一个注意点是数据边界若散斑干涉图存在无信号的黑色边框边界区域应使用掩模处理或保留原始相位图形态不进行填充否则边界填充会把噪声从无信号区域引入有效区域。5. 参数验证与实用技巧让正余弦滤波在工程中真正落地5.1 五招验证滤波结果是否可靠正余弦滤波完成后针对不同的工程目标采用不同验证方式更可靠。下面五个方法覆盖了多数剪切散斑场景的需求解包裹一致性验证。同一张噪声包裹相位图分别用两种不同的解包裹算法解包对比结果。若滤波后两种算法解出的相位差仍超过1.5 rad说明噪声残留无法被解包裹容忍需要增大滤波强度。残差相位图检查。用滤波后的相位去减原始包裹相位再折回区间观察残差图中是否存在条纹或结构。有结构说明滤波强度不足噪声被部分映射成低频假结构。跳变线交叉验证。沿跨跳变线的直线截取滤波前后相位值剖面观察剖面是否仍保持2π折叠行为。若剖面在应跳变的位置出现斜坡则该处滤波核过大。散斑相关指数变化。计算滤波前后散斑相关指数的直方图偏移若直方图右移明显说明有效信号增强。批量比较。在生产线上采集同一物体在不同负载下的多幅剪切散斑图用相同参数滤波观察相位展开结果是否重复性良好。5.2 连续两轮滤波是否值得常见说法是“一次不行就多滤几遍”。实验表明正余弦滤波对相位跳变线的保持能力并非总随次数单调变好。第一轮滤波去除大部分噪声第二轮会把第一轮残余的斑点噪声继续压制但第三轮开始跳变线位置会被反复平均而逐渐模糊。实际生产中用两轮比较稳妥第一轮用sigma1.2的高斯核第二轮用sigma0.8的小核。二轮组合等价于一次中等强度滤波在局部区域的空间分辨率修正比单轮大核保持更多边缘细节。但三轮以上在剪切散斑场景中很少见到正向收益建议做双向消融实验再决定是否用两轮。5.3 与其他预处理步骤顺序的搭配建议正余弦滤波通常放在解包裹之前、质量图计算之后。原因是质量图本身通常由相位梯度构造滤波前计算的质量图能反映原始噪声分布可用于掩模或权重滤波后计算的质量图则反馈的是去噪后的结果用于解包裹阶段更准确。推荐顺序为读取包裹相位图 → 计算初始质量图并生成掩模 → 正余弦滤波 → 更新质量图 → 质量引导解包裹 → 相位展开结果后处理。这套顺序在不同剪切散斑测量系统间基本通用。把正余弦滤波放到解包裹之后属于误用因为在相位展开后的连续相位场上该算法对真实相位梯度的平均会造成不可逆的平滑远不如直接在包裹域处理。5.4 直接嵌入采集管线的脚本将正余弦滤波嵌入处理管线的最终脚本可以参考如下形式def process_shear_speckle(input_path, output_path, kernel_size7, modereflect): # 读取并归一化 data cv2.imread(input_path, cv2.IMREAD_UNCHANGED).astype(np.float32) if data.max() 4: data (data / 65535.0) * 2 * np.pi - np.pi # 正余弦滤波 filtered sin_cos_filter(data, kernel_sizekernel_size, modemode) # 保存位深保持一致 out ((filtered np.pi) / (2 * np.pi) * 65535).astype(np.uint16) cv2.imwrite(output_path, out) process_shear_speckle(speckle_raw.tif, speckle_filtered.tif, kernel_size7)需要根据实际采集软件的输出位深调整65535这个系数。若软件输出为32位浮点相位图则归一化部分可以省略直接使用数值范围判断并做相应裁剪。在实际生产环境中这个函数还应当加入输入尺寸校验、无效值处理与日志记录才能稳定运行在批处理流水线上。正余弦滤波算法本身仅有几十行代码但把它放到适合的位置、调到合适的核大小才能在剪切散斑包裹相位图的后续解包裹与定量分析中真正发挥出价值。本文还有配套的精品资源点击获取