简介基于MATLAB的剪切干涉仪仿真模拟为光学检测与精密测量方向的工程师、科研人员及MATLAB学习者提供可直接复现的数值实验环境该仿真以剪切散斑干涉为核心助力理解物体表面微小不平度、折射率变化及大凹球面缺陷的检测原理。资源压缩包共10个文件主体为8个mexw64编译函数负责光场构建、剪切操作、干涉计算等核心算法辅以1个.m主程序脚本和1张tif示例干涉图整体仅69KB便于下载与快速部署。目前已有2201人学习/下载适合希望借助仿真手段掌握干涉测量原理的光学方向学生与研究者。内容涵盖圆孔光阑、Zernike像差、强度衰减、光束混合等模块通过调用现成函数即可重现两束相干光空间剪切后的干涉图样并进一步分析轴向倾斜、径向弯曲或透明固体内部不均匀性相比从零搭建模型这份资源可大幅缩短理论学习到数值实验的路径。1. 剪切干涉仪的 MATLAB 仿真模拟把波前的“斜率地图”搬到屏幕上在光学检验里横向剪切干涉仪是个让人又爱又恨的结构它不设参考光路把被测波前沿某个方向平移一截再和原波前自己干涉探测器上照样出现条纹。可这些条纹并不是波前形状本身而是波前在剪切方向上的“差分”说白了就是局部斜率。这个反直觉的设计反而让它对高频像差特别敏感所以大口径光学元件检验、自适应光学波前传感里都常见它的身影。这篇文章就用 MATLAB 把这台仪器的仿真模拟完整搭出来从泽尼克波前到干涉条纹再到从条纹反推波前最后给出实际踩坑记录。适合刚入光学测试、正在做 MATLAB 图像处理相关大作业或者想在实验室方案落地前先建个模型的工程师。2. 横向剪切干涉的原理与建模干涉图里藏的是波前差分2.1 用泽尼克多项式描述波前把像差变成 MATLAB 里的矩阵仿真干涉仪的第一步是把被测波前抽象成复数域里的一个二维矩阵。入射光复振幅写成E(x,y) A(x,y) * exp(i * phi(x,y))其中phi 2*pi*W / lambdaW是光程差单位用米或纳米。对于光学系统像差W 最常见的参数化方式是泽尼克多项式展开。泽尼克多项式的口径正交特性让每一项都对应一个特定的像差形态比如离焦、像散、彗差、球差这在光学设计软件里是通用语言。实际做仿真时不需要把全套几百项泽尼克都写出来我一般只保留几个高阶主项用一个简化版 Fringe 泽尼克函数就够了。注意规范化半径 r sqrt(x^2y^2)/R取值 0 到 1孔径外直接置零。代码如下function Z zernikeFringe(j, x, y, R) % 简化版 Fringe Zernike返回单位振幅多项式 % 这里只实现低阶常用项1平移 4离焦 5/6像散 7/8彗差 11球差 % x, y: 网格坐标矩阵R: 孔径半径(像素) theta atan2(y, x); r2 (x.^2 y.^2) ./ R.^2; r sqrt(r2); switch j case 1, Z ones(size(x)); % 平移 case 4, Z 2*r2 - 1; % 离焦 case 5, Z r2 .* cos(2*theta); % 0度方向像散 case 6, Z r2 .* sin(2*theta); % 45度方向像散 case 7, Z (3*r.^3 - 2*r) .* cos(theta); % x方向彗差 case 8, Z (3*r.^3 - 2*r) .* sin(theta); % y方向彗差 case 11, Z 6*r2.^2 - 6*r2 1; % 球差 otherwise, error(未实现的泽尼克项请补表); end end逻辑说明输入网格坐标矩阵和孔径半径 R返回对应项的单位振幅泽尼克面形。代码里用r2代替逐点开方在 512×512 网格上能省掉不少计算量。参数注意两点一是 R 必须小于等于网格边长一半否则孔径会切到网格外二是泽尼克项的归一化约定各家不一这里用的是 Fringe 约定系数含义是该项在孔径边缘造成的最大光程差不是 RMS 值。如果你从光学设计软件导出系数先确认它的归一化口径这一步错了后面所有条纹都会跟着错。2.2 横向剪切的干涉公式平移一份波前和原波前“自己和自己比”横向剪切干涉仪的关键操作只有一个把波前沿着某个方向移动 s 个像素然后和原波前进行复振幅叠加。数学上写成Es(x,y) E(xs, y)叠加后的干涉强度为I |E Es|^2。假设振幅均匀为 1展开得到I 2 2*cos( phi(xs,y) - phi(x,y) )这个相位差就是剪切干涉里最核心的物理量。当剪切量 s 远小于波前变化尺度时phi(xs) - phi(x) ≈ s * d(phi)/dx所以干涉条纹反映的是波前沿剪切方向的梯度而不是波前本身。这就是我开头说的“斜率地图”。反过来说如果我们在波前里塞进一个常数倾斜剪切后这个倾斜会变成常数值的相位差表现为整个条纹图里叠加了均匀的载波频率这个特性在后面单帧条纹提取相位时会用到。MATLAB 里验证这个公式只需要三行deltaPhi phi(:, 1s:end) - phi(:, 1:end-s); % 剪切相位差 I 2 2*cos(deltaPhi); % 干涉强度振幅为1逻辑说明s是剪切量单位是像素。第一行做的是把 phi 右移 s 列后和原场相减得到相位差分图第二行按双光束干涉强度公式生成条纹。参数上要注意剪切之后两个场只在公共区域有干涉意义deltaPhi的尺寸已经比原波前小了 s 列后面做任何处理都要沿用这个裁减后的坐标系不要和原波前的坐标混用。2.3 剪切量的三档选择2 像素、8 像素还是 32 像素剪切量 s 是仿真里第一个要拍的参数它直接决定条纹密度和重建效果。s 太小相位差接近零条纹稀疏到几乎看不出形态对噪声极敏感s 太大斜率被放大条纹密到超过探测器奈奎斯特频率直接混叠。我习惯在仿真里分三档测试剪切量 s像素条纹形态梯度近似程度主要风险2条纹极稀疏近零级观察接近真实导数但信噪比差噪声主导容易重建出一片平8中等密度便于目视和提取线性近似好常用中间档需确认条纹周期大于 4 像素32条纹密边界处容易混叠差分近似变粗趋近“相干微分”高频区折叠必须加密网格参数关联上条纹频率大约等于“波前斜率乘以剪切量除以波长”所以波前本身空间频率越高s 就要取得越小。一个相对稳妥的做法是先按最大波前斜率粗算一次让条纹周期不小于 4 像素再把 s 折半复测。后面第 4 章我会专门讲混叠怎么排查。需要提醒的是s 不必是偶数但如果是单像素奇数值裁剪边界会出现半像素非对称对后续积分重建影响很小可忽略想让边界干净一些就选偶数。3. 用 MATLAB 跑通剪切干涉仿真从波前生成到波前重建一条龙3.1 参数区与波前生成波长、孔径、剪切量写在一个脚本顶部仿真的可复现性依赖参数的集中管理。我会把波长、网格数、孔径半径、剪切量、泽尼克系数全部放在脚本最前面的参数区并给每一项写注释。这样改参数跑对比时不至于在几十行代码里到处翻。% 剪切干涉仪仿真参数区 clear; close all; rng(11); N 512; % 网格边长像素 R 200; % 孔径半径像素 lambda 632.8e-9; % 波长氦氖激光 632.8nm s 8; % 剪切量像素 f0 0.02; % 载波频率cycles/pixel仅用于单帧提取 % 构造网格 [x, y] meshgrid(-N/2:N/2-1, -N/2:N/2-1); rho2 (x.^2 y.^2) / R^2; apertureMask rho2 1; % 圆形孔径掩膜 % 波前光程差 W(m)由离焦 球差组成 coef_defoc 0.5; % 离焦系数单位 lambda coef_spher 1.2; % 球差系数单位 lambda Wl coef_defoc * (2*rho2 - 1) .* lambda ... coef_spher * (6*rho2.^2 - 6*rho2 1) .* lambda; Wl(~apertureMask) 0; % 孔径外光程差置零 % 复振幅加入载波倾斜为后续条纹提取留余量 phi 2*pi .* Wl ./ lambda 2*pi .* f0 .* x; E exp(1i .* phi);逻辑说明这段代码先造一个 512×512 的坐标网格用rho2 1生成圆形孔径掩膜。波前光程差由离焦和球差两项叠加系数以波长为单位所以乘上 lambda 变成米制量纲。最后把光程差转成相位 phi再加上一个 x 方向的线性载波项2*pi*f0*x这个载波在剪切之后会变成固定频率的条纹给第 3 节的 FFT 相位提取提供“单帧即可解调”的条件。参数说明coef_defoc 0.5表示离焦项贡献 0.5 个波长的波前差coef_spher 1.2表示球差贡献 1.2 个波长。球差系数我故意取得比离焦大这样最终波前有一个明显的中心区域便于观察剪切条纹的非均匀分布。如果只想跑通流程coef_spher可以降到 0.5条纹会平缓很多。这个脚本不依赖任何工具箱基础 MATLAB 就能跑不需要找代跑程序之类的外援。3.2 计算剪切干涉图索引切片和 circshift 的边界差异生成干涉图的核心是“平移一份波前再叠加”。很多初学者会直接用circshift(E, [0 s])这其实是坑。circshift是循环移位右移后左边界补过来的数据是原来最右边的数据这在光学上毫无意义等于在剪切场里塞了一圈假信息会在边界形成几条非常亮的人造条纹。正确做法是用索引切片保留公共区域或者干脆把无效区抹成 0。下面给出 x 方向的剪切干涉图计算% x 方向剪切把 E 向右移 s 像素去掉回绕 Ex zeros(N, N); Ex(:, 1:N-s) E(:, 1s:N); % 剪切副本右侧补零 % 有效区域掩膜去掉交界处 s 列 validx zeros(N, N); validx(:, 1:N-s) 1; % 原波前同样裁剪只保留公共区域 Emask E .* validx; Ix abs(Emask Ex).^2 .* validx; % 干涉强度零填充区为0 % y 方向剪切同理 Ey zeros(N, N); Ey(1:N-s, :) E(1s:N, :); validy zeros(N, N); validy(1:N-s, :) 1; Ix abs(E .* validy Ey).^2 .* validy;逻辑说明Ex(:, 1:N-s) E(:, 1s:N)把原始场第 s1 列到最后一列搬到了副本的第 1 列到 N-s 列相当于整体向左平移。注意我这里没有对原场做循环位移而是把副本和原场的公共区求干涉最后用validx把无效区域乘成 0。这样既避免了假信息又保证了Ix和原始网格同尺寸后面频域处理不用再对坐标做一次偏移校正。参数说明s8 时前 8 列和后 8 列的掩膜值为 0也就是说最终有效干涉区是中间的 N-s 列。这条边界让干涉图看起来比孔径小了一圈但波前直径 400 像素远大于 8 像素不影响视觉判断。如果做定量重建重建结果也要使用同样的边界条件避免积分时引入外部伪影。y 方向剪切完全对称只是维度顺序从列变成行。3.3 相位解调单帧条纹图用 FFT 带通滤波提取包裹相位有了干涉图Ix需要把deltaPhi从条纹中解出来。工业上最常用的单帧方法是傅里叶变换法它利用载波把条纹信息搬到频谱里偏离原点的一个峰上用带通滤波器把这个峰抠出来再反变换取辐角。这个过程在光学测量里叫空间载波相移法MATLAB 里实现并不复杂。% FFT 带通相位提取输入干涉强度图 I fftI fftshift(fft2(Ix)); [fx, fy] meshgrid((-N/2:N/2-1)/N, (-N/2:N/2-1)/N); fc f0 * s; % 剪切后载波频率 原载波频率 * 剪切量 bw 0.04; % 带通高斯宽度单位 cycles/pixel % 带通滤波器中心在 (fc, 0) H exp(-((fx - fc).^2 fy.^2) ./ (2*bw^2)); analytic ifft2(ifftshift(H .* fftI)); wrapped_x angle(analytic); % 包裹相位范围 [-pi, pi)逻辑说明fftshift(fft2(Ix))得到频谱频谱里原点和正负载波频率处各有一个峰。带通滤波器H是中心在 (fc, 0) 的高斯窗只保留正载波峰抑制背景和负频峰。反变换得到的analytic是一个复解析信号取辐角就是包裹在 [-pi, pi) 里的相位差deltaPhi。之所以能这么干是因为条纹强度在数学上可以写成I 2 cos(deltaPhi)而exp(i*deltaPhi)的信息完整地映射到了正频峰里。参数说明fc f0*s 0.02*8 0.16周期/像素这就是条纹载波在频域里的坐标。设计这个参数时至少要满足两个条件fc远大于波前相位梯度的频谱宽度又远小于奈奎斯特频率 0.5。bw 0.04是我常用的起点它决定滤波器的频带宽度太窄会切掉波前的高频细节太宽会把背景峰和负频峰卷进来。遇到被测波前有较大彗差这类空间高频时把bw适当加大到 0.06~0.08再观察条纹边缘有没有“糊掉”。这一步本质上是图像处理里的频域滤波操作和你用 MATLAB 做条纹图片处理时的套路完全一样。解包裹这一步仿真里可以用一维展开加中值对齐来近似。严格说二维解包裹应该用质量引导或最小二乘算法但对仿真验证一维逐行展开并扣除载波已经足够。核心代码% 逐行一维解包裹再减去载波项 d diff(wrapped_x, 1, 2); d d - 2*pi*round((d pi) / (2*pi)); % 把差分值包到 [-pi, pi) ph_unwrap cumsum([wrapped_x(:,1), d], 2); ph_unwrap ph_unwrap - 2*pi*fc .* x; % 减去载波相位 % 得到单位为“米”的波前梯度 dWdx ph_unwrap ./ s ./ (2*pi) .* lambda;逻辑说明第一行对相位差做横向差分第二行用round((dpi)/(2*pi))求出整数个 2π 跳变并扣掉这就是一维解包的原理。cumsum沿行方向积分得到连续相位。随后减去已知载波2*pi*fc*x剩下的就是纯剪切相位差。除以剪切量 s 再乘 lambda/(2π)得到每像素上的波前斜率单位米每像素。这套代码没有用任何工具箱函数全是基础数组运算适合做一个可移植的仿真骨架。3.4 波前重建从两个正交方向的斜率反推波前只从一个方向的剪切只能得到一维斜率要重建二维波前必须有 x 和 y 两个正交剪切方向的干涉图。对 x 方向用Ix得到dWdx对 y 方向用Iy得到dWdy然后用最小二乘积分把它们组合成一个波前。频域积分法很直观空间域里波前梯度等于波前的偏导对应到频域就是乘2πi*u和2πi*v反解一个线性方程即可function Wrec integrateFFT(dWdx, dWdy, lambda) % 频域最小二乘积分重建波前 % dWdx, dWdy: 每像素波前斜率单位 m/像素 [ny, nx] size(dWdx); % 频率坐标单位 cycles/pixel u (-nx/2 : nx/2-1) / nx; v (-ny/2 : ny/2-1) / ny; [U, V] meshgrid(u, v); % 频域微分算子d/dx 对应 2πi*U % 构造并求解 (U.^2V.^2) * F(W) U*F(dWdx) V*F(dWdy) D (2*pi*U).^2 (2*pi*V).^2; D(1,1) 1; % 避免除零 S -(1i*2*pi*U) .* fft2(dWdx) -(1i*2*pi*V) .* fft2(dWdy); Wrec real(ifft2(S ./ D)); Wrec Wrec - mean(Wrec(:)); % 去掉整体平移 end逻辑说明频域里波前梯度dWdx的傅里叶变换等于2πi*U*F(W)dWdy对应2πi*V*F(W)。把两个方向的梯度方程合起来用最小二乘的形式解出F(W)。分母D在零频处置 1既保证不除零又相当于把常数项平移量钳制为 0。最后real(ifft2(...))得到波前减均值是去掉活塞项。这个算法对边界形状没有任何要求圆形孔径也能直接处理比空间域逐行积分干净得多。参数说明注意dWdx和dWdy必须经过有效掩膜裁剪且尺寸一致。如果两个方向的剪切干涉图有效区不同积分前要先用min求交集否则频域积分会把不一致的边界当成突变梯度重建出来会出现一道裂缝。lambda在这里只用来单位统一实际计算时只要dWdx/dWdy量纲一致即可。4. 剪切干涉仿真避坑条纹混叠、方向错位和边界伪影的排查4.1 条纹密度过载导致混叠现象干涉图高频区域出现摩尔纹一样的折叠条纹原本应连续变化的条纹在边缘变成了一圈一圈的假环FFT 提取时该处的相位明显断裂。原因剪切量 s 和局部波前斜率乘积对应的空间频率超过了采样极限。常出现在球差波前的边缘那里相位梯度很大而我把 s 直接从 8 跳到 32 时忘了复核。解决将 s 调小或者把网格 N 翻倍。仿真里最简单的做法是保持 s 不变把 N 从 512 加到 1024孔径半径 R 同步翻倍这样同一切割量对应的物理空间频率不变但采样率提高了。同时检查fc f0*s是否小于 0.4给奈奎斯特留出余量。4.2 重建波前的像散方向颠倒现象明明我在波前里只加了 x 方向的像散项重建结果却显示像散旋转了 90 度或者出现在 y 方向上。原因x 方向剪切和 y 方向剪切得到的梯度在积分时被放反了。常见于把Iy转置后送入integrateFFT或者在载波扣除时把 x 和 y 的载波频率对调了。更多时候是dWdx和dWdy在网格维度上的行列约定不一致。解决在重建前画一个向量场检查用quiver(x(1:8:end,1:8:end), y(1:8:end,1:8:end), dWdx(1:8:end,1:8:end), dWdy(1:8:end,1:8:end))看一眼箭头方向是否为孔径中心向外发散。如果是涡旋状或整体偏转 90 度那基本就是 x/y 梯度互换或转置修正后再积分。4.3 重建波前边缘“翘起”并叠加全场倾斜现象重建出来的波前沿孔径边缘翘得很高中间区域形态正常但整体带一个平面倾斜。原因剪切干涉本身对低频倾斜不敏感但重建积分时边界截断引入了一阶伪差另一个常见来源是载波扣除不准确残留的线性相位混进了梯度。解决先在频域积分后做一次平面拟合Wrec Wrec - polyplane把a*xb*yc的低阶项去掉。更本质的办法是同时用两个不同的剪切量 s1、s2 各重建一版两版结果求差差面型如果仍然翘边说明问题源在掩膜或载波扣除而不是积分算法。4.4 中文字符乱码导致脚本直接报错现象MATLAB 打开 .m 文件中文注释显示为乱码点运行时提示 “Invalid text character” 或直接定位到注释行报错。原因不同版本 MATLAB 对 .m 文件编码的默认值不一致。老版本默认 GBK新版本默认 UTF-8。你用 R2023b 里写的带中文注释脚本拿到别的机器打开编码错位就会触发解析错误。解决统一用英文注释或者把脚本用 UTF-8 编码重新保存。MATLAB 编辑器里“编辑器”选项卡下找到“保存文件编码”选 UTF-8 再保存一次。命令行里执行feature(DefaultCharacterSet,UTF-8)可以临时切换运行环境的字符集但对已有文件的修复还是以重新保存为主。仿真代码本身逻辑不复杂我后来习惯变量名和注释全用英文彻底避开这套坑。4.5 新版 MATLAB 启动闪退或并行池报错现象安装的是 R2026a/R2026b装好打开就闪退或者跑脚本时parpool初始化报错导致没法进入仿真环境。原因新版 MATLAB 在部分显卡驱动和 OpenGL 组合下会闪退并行计算工具箱默认启动并行池时又容易与多核心环境冲突两个问题常一起出现。解决先用matlab -nodesktop -nosplash启动命令行形式能进去就说明图形界面相关配置有问题。进入后在“预设项-常规-图形硬件加速”里改为软件加速。并行池报错的话在代码里显式关闭自动并行parpool(local, 0)或注释掉所有parfor。注意这个坑和环境变量有关不是代码问题别花太多时间在找程序逻辑上。4.6 circshift 回绕制造的边界伪条纹现象干涉图右边界出现一条笔直的高频亮带FFT 提取后这列数据相位异常重建波前在对应位置出现一道纵向裂缝。原因直接用了circshift(E, [0 s])圆周移位把右侧溢出的数据绕到左侧在剪切相位差里产生了一个巨大的假跳变等价于在边界塞了一根高斜率楔形。解决改用本文 3.2 节里的索引切片方式并显式乘以有效区掩膜。如果因为其他原因必须用circshift那在叠加前把E和剪切副本都乘以一个不包含回绕区的掩膜让无效区域归零。这个坑在仿真初版几乎必踩我把它排在避坑清单最后一条是因为它最隐蔽表面看干涉图形态还挺正常只有做定量精度分析时才暴露。5. 把仿真结果验证到可信三步自检与一组交叉验证仿真跑通只是第一步真正能投入到方案验证需要确认重建波前和输入波前高度一致。我最常做的第一个自检是解析梯度对比把输入的离焦和球差函数求解析偏导和dWdx直接做点对点比较。离焦项W c*(2*rho2-1)对 x 的偏导是4*c*x/R^2球差项偏导也不难写把这个解析梯度作为基准算出重建梯度的 RMS 误差。在 N512、s8、无噪声条件下RMS 误差应该在 λ/100 量级如果差到 λ/10一定是哪个环节出了系统性问题。第二个自检是剪切量无关性测试。剪切量是算法里的可调参数不是物理系统的本征量因此用 s4 和 s16 分别重建同一个输入波前两次结果的差面型应该接近零。如果两版波前差异超过 λ/50说明算法对剪切量的依赖过大通常是载波参数或滤波器带宽没配对。这个测试只需要改一个数字跑两组脚本成本极低我每次写完新的仿真流程必跑这一条。第三个自检是相位提取算法的等效性检验。我用 FFT 法跑完一组再用三步相移法在仿真里生成三帧相移条纹重新提取相位两套流程的重建结果应该一致到小数点后两位。其实仿真里的“相位提取”本质上是在验证你自己的理解不是验证算法优劣所以只要关注结果差面型即可。三步相移法在仿真里实现很直接给波前基底分别加 0、2π/3、4π/3 的常量相位漂移生成三帧干涉图用反正切公式就能恢复出包裹相位。整个过程不用真实验光路跑起来非常快可以用来反查 FFT 法里载波频率设置是否有偏差。最后补一个我个人很依赖的习惯把所有重建结果保存成标准格式包括输入波前、重建波前、残差面型和 RMS/PV 数值。跑过几次参数扫描之后你会发现很多看起来“玄学”的误差波动其实都来自某一组固定参数组合而不是随机噪声。保留残差面型比只存数值有用得多因为残差图能直接告诉你误差集中在孔径边缘还是中心是载波泄漏还是剪切边界问题。这套仿真骨架同样可以扩展到底面形误差检测、大口径拼接检验方案预演的思路里核心就是把“剪切干涉”这个物理过程还原成矩阵运算把每一处边界条件都管住剩下的交给 MATLAB 的数组运算就行。希望帮到你。本文还有配套的精品资源点击获取