MATLAB实现达曼光栅设计:从GS相位搜索到等强度点阵仿真

📅 发布时间:2026/10/13 20:00:30
MATLAB实现达曼光栅设计:从GS相位搜索到等强度点阵仿真
简介在光学实验与激光整形中达曼光栅可将单束光转换为等间距点阵而MATLAB模拟是观察其干涉效应的高效途径。这份压缩包面向光学专业学生、科研人员及MATLAB数值仿真学习者包含一个完整的M脚本通过定义光栅周期、透镜尺寸与相位延迟等参数结合二维傅里叶变换提取衍射峰并以图像形式直观呈现点阵分布。压缩包仅719字节内含1个M文件代码量虽小却覆盖了从参数设置、相位计算到频域分析和可视化的核心流程适合用于课程设计或快速验证光学理论。该资源已有2324人学习作者为beau_lily程序结构清晰便于用户逐步调试并理解达曼光栅的设计逻辑通过运行和修改参数还能进一步探索光栅周期、相位延迟对点阵密度与干涉效果的影响是衔接光学原理与工程仿真的实用工具。1. 达曼光栅生成点阵把一束激光拆成等强度网格的标准做法如果你做激光加工、光镊阵列或者并行光互连迟早会遇到一个需求一束高斯光进来出射端要变成 3×3、7×7 甚至 11×11 的等强度点阵。用透镜阵列太重用空间光调制器太贵而达曼光栅恰好是这类分束场景里最常见、也最成熟的光学方案。所谓达曼光栅本质上是一层二元相位片只含 0 和 π 两种相位台阶设计得当就能让入射光均匀分配到指定的一组衍射级上远场看到的就是整齐的点阵。这篇笔记要讲的就是怎么用 MATLAB 把“设计达曼光栅相位轮廓”这件事跑通从相位搜索到点阵仿真再到底图导出和常见排错全程可以照着复现。2. 相位模型与傅里叶展开达曼光栅“复制光斑”的数学基础2.1 二元相位结构为什么 0/π 台阶能分束达曼光栅的透过率函数写出来并不复杂t(x) exp(j·φ(x))其中 φ(x) 只取两个值0 和 π。当光穿过一个周期性的相位结构远场光强分布由透过率函数的傅里叶展开决定每一项展开系数对应一个衍射级光强就是系数模的平方。这样设计达曼光栅的问题就转化成在一个周期内安排若干段 0/π 相位台阶让目标衍射级的傅里叶系数振幅尽可能相等同时让非目标级次的能量尽可能少。实际操作中相位 0/π 对应复平面上的 1 和 -1所以二值相位结构等价于一个实数振幅光栅这比连续相位要更容易加工和复制。你拿扫描电镜看一片达曼光栅看到的就是一条条矩形沟槽槽深对应 π 相位差槽宽和位置则决定了能量分配。这也是达曼光栅比全息分束器更受工业欢迎的原因二元轮廓对工艺误差的容忍度比连续浮雕高很多。当然并非随便排列 0/π 台阶就能得到均匀点阵。相邻两次相位翻转之间的宽度、翻转点的位置都直接影响各衍射级的能量分配。要找到一组最优翻转点就是设计程序的核心任务也就是标题里那个 MATLAB 程序真正在做的事。2.2 二维点阵如何合成从一维序列到矩阵外积一维达曼光栅只能把光分散到一条线上的若干衍射级。要得到二维点阵常见做法是把两个一维达曼光栅的方向垂直叠加这在数学上表现为相位矩阵的外积相加。假设一维相位序列为一个列向量 φ(x)共 N 点只含 0 和 π。二维相位矩阵按下式构造φ2(i,j) φ(i) φ(j)由于 exp(jπ·整数) 只有 1/-1 两种取值二维透过率矩阵实际上就是一维透过率向量的外积T2 t(x) ⊗ t(y)换句话说二维达曼光栅的衍射级能量谱是一维级次谱的二维卷积展开。如果你设计了 9 个一维衍射级二维合成后就是 9×9 81 个点。这个外积性质非常有用设计二维光栅时先搜一维相位序列再直接合成比直接在二维平面上做优化快一个量级。仿真验证时二维远场强度由二维傅里叶变换得到离散实现就是在 MATLAB 里对 T2 做 fft2。这一步计算量很小一个 128×128 的矩阵毫秒级就能跑完适合反复调参数。2.3 两个核心指标均匀性与衍射效率达曼光栅的好坏业内主要看两个指标。第一个是均匀性误差定义为点阵中各点光强的最大偏差与平均光强的比值。理想情况下所有点等亮均匀性误差为 0实际设计做到 3% 以内已经算优秀做到 1% 以下则需要非常精细的初值和足够的迭代次数。第二个指标是衍射效率指的是落在目标级次上的能量占总透射能量的比例。达曼光栅的衍射效率理论上可以做到 80% 以上但均匀性和效率往往此消彼长过度压均匀性会牺牲效率过度追求效率又会让个别点明显偏暗。实际工程中我会建议先固定目标效率底线再最小化均匀性误差而不是两个目标同时优化。GS 算法Gerchberg-Saxton是搜索这类二元相位轮廓最常用的起点。它的思路很直白在空间域约束相位为 0/π在频域约束目标级次振幅相等反复迭代直到收敛。下一章就直接进入 MATLAB 实现。3. 在 MATLAB 里跑通最小程序GS 迭代求解一维相位并合成二维点阵3.1 一维 GS 迭代用投影矩阵直接算各级次振幅一套可复现的最小程序我习惯把 DFT 写成投影矩阵而不是直接调 fft。原因很简单fft2 的频率坐标受采样点数限制级次未必正好落在整数频点上投影矩阵则直接针对我们关心的那些衍射级次计算傅里叶系数物理含义更干净也方便后续算均匀性。以下是一维相位搜索代码目标为 9 个等强度级次% 一维达曼光栅 GS 相位搜索 N 128; % 一个周期内的采样点数 orders -4:4; % 目标级次 -4 到 4共 9 级 Mlen length(orders); % 目标级次数 % 频域投影矩阵A(m,n) exp(-j*2*pi*orders(m)*n/N) / sqrt(N) A exp(-1j * 2 * pi * (orders. * (0:N-1)) / N) / sqrt(N); target ones(Mlen, 1); % 目标级次振幅全部为 1 phase 2*pi*rand(N, 1); % 随机初相位 for iter 1:300 % 前向当前空间相位对应的各衍射级复振幅 E A * exp(1j * phase); % 频域约束振幅强制为1保留相位 E target .* exp(1j * angle(E)); % 逆向由目标频谱反推空间相位 phase angle(A * E); % 去掉整体相位偏置避免数值发散 phase phase - mean(phase); end % 二值化到 0 / pi phase_binary pi * (phase 0); % 用二值化相位重新计算各级次强度 E_final A * exp(1j * phase_binary); I_1d abs(E_final).^2; disp(I_1d.);这段代码里A 是 9×128 的投影矩阵做的事情就是把 N 点空间信号投影到 9 个衍射级上。每次迭代分三步先算当前相位的频域分布再把目标级次的振幅改成 1最后把修正后的频谱反推回空间相位。循环 300 次后把连续相位压到 0/π 两个值上得到可直接用于加工的一维相位序列。两个参数值得注意。N128 是每周期采样数它决定相位台阶的位置精度128 对 9 级光栅足够迭代次数 300 属于保守设置多数情况下 100 次就能看到级次强度趋向一致再多也不会明显改善。随机初值会影响最终相位序列长什么样但不影响点阵质量这是 GS 迭代本身的特点。3.2 二维合成与点阵显示外积拼接加 fft2拿到一维相位序列后按外积规则生成二维相位矩阵然后用 fft2 观察点阵。为了让最终画面更接近真实远场需要把单个周期平铺若干倍再做变换这样频谱上出现的是锐利的峰值而不是糊成一片的谱线。% 二维达曼光栅合成与点阵显示 t1 exp(1j * phase_binary); % 一维透过率向量 T2 t1 * t1.; % 外积得到二维透过率矩阵 % 用投影矩阵直接算 9x9 个级次的复振幅精确计算 E2 A * T2 * A.; % 9x9 复数矩阵 I2 abs(E2).^2; % 各点光强 % 显示用把基元周期平铺 20x20模拟真实周期结构 G 20; T2_tiled repmat(T2, G, G); spec fftshift(fft2(T2_tiled)); power abs(spec).^2; % 截取中心区域显示点阵 imagesc(power(600:1000, 600:1000)); axis image; colormap hot; title(达曼光栅模拟远场点阵);投影矩阵计算出来的 E2 是 9×9 复数矩阵直接取模平方就得到了 81 个点光强的相对值这是评估均匀性最可靠的数据源。而repmat(T2, G, G)是为了出图好看它相当于在空间上重复摆放基元周期模拟真实的周期性光栅结构。fft2 之后用 fftshift 把零频移到中心才能看到对称的点阵。T2_tiled 是 2560×2560 的矩阵在普通笔记本上做一次 fft2 大概一两秒完全在可接受范围。如果你只需要数据不需要图可以跳过平铺部分直接信任 E2 的结果。3.3 均匀性和效率的定量计算肉眼看到 9×9 亮点均匀还不够必须给出数值。均匀性和衍射效率的统一计算方法如下% 定量评估 I_vec I2(:); % 81 个点的光强 uniformity (max(I_vec) - min(I_vec)) / mean(I_vec); fprintf(均匀性误差%.2f%%\n, uniformity * 100); % 衍射效率目标级次能量占比 powTotal sum(abs(T2(:)).^2) / N^2; % 入射总能量归一化 eff sum(I_vec) / powTotal; fprintf(衍射效率%.2f%%\n, eff * 100);这里 powTotal 用了帕塞瓦尔定理空间域平均能量等于频域总功率。因为我们用的是归一化投影矩阵所以 E2 的功率总和可以直接和 T2 的平均功率做比。实际跑下来一个 9×9 的达曼光栅程序效率通常在 70% 上下波动均匀性在理想迭代条件下能压到 2% 以内。一个容易踩的误区是直接用max(I2(:)) / min(I2(:))当均匀性指标这个比值在理想均匀时是 1但实际评估时对个别暗点过于敏感不如最大偏差相对均值直观。上面代码用 (max-min)/mean既有明确物理含义也和加工容忍度直接对应。4. 像素数、级次与波长的取舍仿真前必调的 3 组参数4.1 采样点数 N 与最大偏转角Nyquist 条件不能破一个周期内的采样点 N 直接决定了可达的最大衍射角。设空间采样步长为 Δ衍射级次 m 的偏转角满足sinθ ≈ m·λ / (N·Δ)要保证最高级次能正常衍射必须满足 m·λ/(N·Δ) 1。N 太小高级次会进入渐逝波区域仿真结果直接失真表现为点阵边缘的级次明显偏暗。所以 N 的选择不是越高越好而是要匹配目标级次数和实际工艺最小线宽。常见做法是取 N ≥ 4M比如 9 级光栅 N128已经远大于 36。N 偏大时相位台阶的宽度仍在整数个像素上对 fft 计算没有本质影响但投影矩阵 A 变成 N×M 的稠密矩阵迭代耗时随 N 线性增长。N 128 到 256 之间300 次迭代都只是一两秒的事不必刻意省。4.2 目标级次 M 取奇数还是偶数奇数是默认选择达曼光栅的目标级次几乎总是奇数因为需要保持衍射级对称分布这样相位序列本身也相对容易搜索。下面是不同级次数的典型参数对比目标级次 M衍射级分布一维点数二维点阵数单点相对能量典型应用3-1,0,139高少光束并行加工5-2,-1,0,1,2525中光镊阵列、激光打标7-3,-2,-1,0,1,2,3749中低多通道并行探测9-4,-3,-2,-1,0,1,2,3,4981低大面积面阵照明11-5,-4,-3,-2,-1,0,1,2,3,4,511121很低高密度阵列M 越大每个衍射级分到的能量越小。举例来说3×3 点阵每个点能拿到约 20% 的入射能量而 11×11 点阵每个点只有约 6%这是在选择级次数时必须接受的物理现实。程序上只需修改 orders 变量的范围其他代码不用动。4.3 波长、像素尺寸与相位深度的换算GS 程序输出的相位是 0/π 的抽象值真正加工时要把 π 相位差换算成刻蚀深度。换算公式是h λ / (2·(n - 1))其中 λ 是工作波长n 是基底材料折射率。以石英在 532 nm 波长为例n≈1.46h≈532/(2×0.46)≈578 nm。这个深度只和波长、折射率有关和光栅周期无关所以同一套 GS 程序换波长时只需要改 h不需要重新搜索相位。像素尺寸 Δ 则决定了光栅的物理周期 d N·Δ也就决定了相邻衍射级的角间距Δθ ≈ λ / d如果你需要在某个工作距离上得到特定间隔的点阵先确定 Δθ反推 d再由 d 和 N 确定 Δ。这是从仿真走向光路调试时最实际的一步很多新手在 MATLAB 里看图形很漂亮一到光学平台就对不上点阵间距原因就是漏了这一层换算。4.4 迭代次数与初值快收敛的经验区间GS 迭代对初值不挑剔但收敛速度差异很大。我一般用三种初值做对比随机二进制相位、随机连续相位、以及把 0/π 按 Fibonacci 序列粗排的硬限幅初值。随机连续相位收敛最平滑不容易陷入局部振荡随机二进制初值如果配合二值化约束收敛速度快但容易卡在次优解。迭代次数的经验区间是 100 到 500。100 次之后均匀性曲线基本不再下降这时可以提前跳出循环。建议在循环里每 20 次打印一次均匀性观察曲线是否还有下降趋势。如果连续 3 次打印结果不变直接 break不必死等。5. 达曼光栅程序落地避坑5 个常见问题与排查5.1 二值化之后均匀性突然变差现象GS 迭代过程中连续相位算出来的级次强度非常均匀误差不到 1%但最后强制二值化到 0/π 后均匀性跳到了 8% 以上。原因连续相位有更丰富的调节自由度能够把能量微调到每个级次上二值化相当于把空间相位锁死到两个离散值高频分量被重新分配各衍射级能量自然出现偏差。这不是程序 bug而是二元光学设计固有的代价。解决不要用连续相位的评估结果作为最终指标始终用二值化后的相位重新计算。如果二值化后均匀性超标有两条路一是把二值化约束提前到 GS 循环内部每次迭代后直接把相位压到 0/π 再继续二是改用模拟退火或遗传算法做翻转点搜索它们能直接以二值相位为自由度结果往往比 GS 二值化后更好。多数情况下我会先试第一种改动小跑得快。5.2 仿真点阵出现十字亮线主点周围糊成一片现象用 fft2 直接变换 T2 后点阵四周有一整条亮十字线目标点之间还有明显的不该出现的亮斑。原因直接用单个周期做 fft2相当于把一个非周期截断的边界当成了硬孔径频谱里混入了矩形窗的 sinc 拖尾。T2 矩阵边缘与矩阵中心对应的是同一周期内的不同位置但 DFT 默认信号是周期延拓的平铺倍数小时边界效应就会呈现为十字形条纹。解决出图前用repmat把基元周期平铺 10 到 20 倍再做 fft2。平铺相当于让信号更接近真实周期结构谱峰更锐利拖尾被摊薄。如果图像依然有微弱十字线考虑加 Hanning 窗做切趾但要知道这会改变光强比例不适合用于指标计算。5.3 fftshift 忘了用点阵被劈成四块现象点阵图像看起来不是中心对称的方形网格而是四个角落各有一小簇亮点中心是黑的。原因MATLAB 的 fft2 输出的零频在矩阵左上角不做 fftshift 就直接显示零级跑到左上角其他各级按周期性折叠到四角看起来就是四块。解决显示功率谱之前先fftshift(fft2(...))。这个错误几乎每个初学衍射仿真的人都会遇到一次不算丢人但会影响后续对点阵坐标的判断。在验证代码里把 fftshift 写死在显示段后续改代码就不容易漏。5.4 偶数级次设计时点阵中心永远是暗的现象尝试做 2×2 或 4×4 的等分点阵把 orders 设置成 -2,-1,1,2 这类不含 0 的对称级次搜出来的相位序列在仿真中总是对应不上中心区域出现暗点。原因偶数级次要求抑制零级但 GS 算法的频域约束里如果没有明确把零级振幅设为 0迭代会自然地把能量留在零级附近。即使约束了振幅为 0二值化后零级往往还是压制不干净物理上这与二元光栅的结构对称性有关。解决尽量避免偶数级次的达曼光栅设计。需要 2×2 或 4×4 分束时可以用奇数级次的组合比如设计 5×5 再加掩膜挡住多余点或者用两个一维奇数光栅差频实现。如果你坚持要偶数级次需要把零级约束权重提高并接受均匀性劣化的代价。工程上我几乎没见过偶数级次达曼光栅做得好用的常规方案都会绕开。5.5 换了 MATLAB 版本或加了并行池结果不一样现象同一套代码在旧版 MATLAB 跑均匀性 1.5%在新版跑变成 2.8%或者开了并行池之后某几步循环结果出现微小差异。原因GS 迭代依赖随机数初值rand函数的随机流在不同版本、不同启动状态下生成的序列不同导致最终相位序列不同。二值化后若干翻转点边界处于临界位置时微小数值扰动会让相位多跳一格强度分布随之偏移。解决固定初值。在代码开头加rng(2024)这类随机种子设置保证每次运行得到同一结果。如果需要在不同电脑间对比把最终相位序列存成 .mat 或 .csv 文件作为共享基准而不是每次现跑现比。这个习惯能省掉大量“我这边明明调好了”的沟通成本。6. 从仿真到加工相位轮廓导出灰度图与版图检查的实用细节6.1 相位矩阵转灰度图注意周期闭合和线宽检查达曼光栅最终要变成光刻版图或直写文件最常见的形式是 8-bit 灰度图。0 相位对应黑色π 相位对应白色灰度图按 256 级输出。导出前要做两件事确认相位矩阵左右、上下边缘是闭合的因为光栅是周期结构边缘处相位不连续会导致实际加工出来的器件边界处出现多余衍射另一个是检查最小线宽。% 导出灰度图并检查周期闭合 % 强制首尾行/列一致保证周期闭合 phase_closed phase_binary; phase_closed(end, :) phase_closed(1, :); phase_closed(:, end) phase_closed(:, 1); % 归一化到 0-255 img uint8(255 * (phase_closed 0)); imwrite(img, dammann_9x9.png); % 最小台阶宽度统计相位翻转点之间的间距 flipIdx find(diff(phase_binary) ~ 0); minLine min(diff(flipIdx)); fprintf(最小台阶宽度%d 个像素\n, minLine);最小线宽直接决定能不能被光刻工艺做出来。假设像素尺寸是 1 μm最小台阶宽 3 个像素就是 3 μm一般紫外光刻能胜任如果最小线宽只有 1 个像素就需要换更高分辨率的直写设备成本立刻上一个台阶。所以看到 minLine 太小时优先考虑增大 N 或改用更平滑的相位序列而不是硬着头皮送加工。6.2 相位深度公差加工前先做一个容差仿真相位深度误差是达曼光栅实际器件性能波动的主要来源之一。刻蚀深度偏了 5%各衍射级强度比例就会整体变化。稳妥的做法是在交付版图前把相位偏移量设置为 ±5%、±10% 各跑一组完整仿真看均匀性是否会劣化到不可接受。我在实际项目中遇到过相位深度偏了 8%9×9 点阵边缘点明显变暗的情况。原因是刻蚀速率受图形密度影响沟槽宽的地方刻得深窄的地方刻得浅。现在多数光刻代工厂能提供工艺补偿但你至少要在设计阶段给出允许公差范围把补偿量写在技术协议里。6.3 拿到实物后的验证顺序器件到手先别急着上万能支架调光路。用功率计在焦平面上逐点测 9×9 个点的光强计算实测均匀性和仿真值对比。如果偏差在 3% 以内说明设计、加工、对光三环节都没问题如果偏差明显先检查入射光斑是否完全覆盖光栅有效口径再来回微调入射角。很多时候器件本身没问题是入射光没对准或者光斑大小和光栅通光口径不匹配。我自己养成的习惯是任何一次设计都保存三份东西随机种子、相位矩阵原始 .mat 文件、二值化灰度图。这样无论是换电脑复现还是加工厂反馈某个细节问题都能快速定位是设计问题还是流程问题。这个“后悔药”成本极低但能避免很多重复工作。希望帮到你。本文还有配套的精品资源点击获取