★ FEATURED ARTICLE 简介《光学工程MATLAB开发超表面透镜逆向设计优化算法》是一份面向光学工程与计算成像领域的高阶技术笔记主要服务于需要掌握超表面器件设计与优化算法的科研人员、工程师及研究生。文档从光学工程与超表面透镜的基本概念入手系统讲解逆向设计的数学基础涵盖麦克斯韦方程组、波动光学理论与相位调控机制并重点展开遗传算法、粒子群算法和模拟退火算法在MATLAB中的目标函数构建、参数设置、迭代收敛与结果可视化等完整实现步骤配有高分辨率成像、多焦点和宽带透镜等实际案例的性能评估与对比。资源包仅含1个PDF文件大小约4.11MB支持目录章节跳转与阅读器大纲快速定位文字、图表、函数与目录均显示正常内容完整、条理清晰。该资料已有103人学习浏览适合希望从理论到代码系统掌握超表面透镜逆向设计全流程的读者参考使用。1. 超表面透镜逆向设计为什么“正向试凑”在亚波长结构上行不通超表面透镜由亚波长纳米柱阵列构成通过局部相位延迟整形波前。很多人一开始按传统透镜思路走仿真几个柱径看焦点改柱径再仿真一两天下来结果往往只是局部改善NA稍微拉高就失控。问题在于结构自由度太多互相耦合焦点质量对每个单元都不敏感人工试错没有全局视角。逆向设计正好反着来把“想要的光场分布”写进目标函数用优化算法自动迭代每个单元该长成什么样解决的是从目标光学响应反推结构参数的反问题。在光学工程里这套方法覆盖大数值孔径超表面聚焦、多波长消色差、偏振复用以及特殊焦斑形态的波前整形。MATLAB 在这个流程里的位置比较特殊它不是电磁仿真主力FDTD/RCWA 通常在外部却是把“目标函数前向传播伴随梯度优化循环后处理”整合成闭环最顺手的环境调试和参数实验成本低。适合读这篇的人是正在用 MATLAB 做微纳光学设计的光学工程师以及刚接手超表面项目、想拿到第一份可跑通代码的从业者。2. 正向模型与设计目标逆向算法落地前的地基2.1 从焦距到目标相位先知道透镜该长什么样逆向设计的第一件事是把正向过程写得足够简单可靠。对超表面透镜而言正向过程的核心是相位分布到焦面光场的传播。目标焦距、设计波长、透镜半径确定后理想相位分布只有一个解析解双曲相位。设径向坐标 r焦距 f波长 λ相位延迟为φideal(r) −(2π/λ)·(√(r²f²) − f)这个式子有两点要强调相位不是 r² 的简单二次型符号约定必须和后续透射系数公式一致。我一般把单元透射系数写成 t(r) A(r)·exp(jφ(r))φ 是相对自由空间传播多出来的相位延迟所以 φideal 是负值。很多代码第一步符号设反后面梯度符号全反损失下降但优化出一个“反相位透镜”。用代码封装时直接算精确式比用近轴近似更稳。半径 10 μm、焦距 20 μm、λ633 nm 时边缘相位约 −23.4 rad近轴近似 −πr²/(λf) 得 −24.8 rad差 1.4 rad已超过 π/2对焦斑质量的严格要求来说不可忽略。function phi target_phase(r, lambda, f) % 计算超表面透镜的理想双曲相位分布 % r : 离轴半径, 单位 um, 可为标量/向量/meshgrid矩阵 % lambda: 设计波长, um % f : 设计焦距, um % 返回值 : 相位延迟, rad, 负值 phi - (2*pi/lambda) .* (sqrt(r.^2 f^2) - f); end这个函数有两处用途作为优化算法的初始解避免随机初始相位造成的不必要局部最优在评估阶段计算理想波前差。注意 r 可以直接传 meshgrid 结果MATLAB 的数组运算比 for 循环快一个数量级。理想相位定完后还要定义口径口径外的相位不影响聚焦却会在传播中产生边缘衍射优化结果容易在口径边缘堆积异常相位应加一个圆形掩模把设计变量限制在口径内并在边缘留 2~3 个网格过渡带。2.2 把前向传播封装成MATLAB函数角谱法的关键实现透镜面复振幅经过相位调制后要传播到焦面计算光强。前向传播算子我固定用角谱法而不是 Fresnel 积分。原因有两个角谱法对传播距离没有近轴限制它是酉变换伴随算子可以直接复用同一段传播代码只是方向反过来。这一点对后面伴随梯度实现很重要。function E_out prop_asm(E_in, z, lambda, dx) % 角谱法衍射传播算子 % E_in : 入射复振幅, Ny x Nx % z : 传播距离, um; z0 表示反向传播伴随用 % lambda: 波长, um % dx : 网格间距, um [Ny, Nx] size(E_in); fx (-Nx/2 : Nx/2-1) / (Nx*dx); % 空间频率, 1/um fy (-Ny/2 : Ny/2-1) / (Ny*dx); [FX, FY] meshgrid(fx, fy); k 2*pi/lambda; kz sqrt(k^2 - (2*pi*FX).^2 - (2*pi*FY).^2); kz(imag(kz) ~ 0) 0; % 倏逝波置零 H exp(1j * kz * z); % 频域传播相位 E_out ifft2(fft2(E_in) .* ifftshift(H)); end这段代码有三处易错点。第一频率网格 fx 范围是 [−1/(2dx), 1/(2dx)]单位 1/um必须与 dx 的物理单位一致把 dx 的尺寸定义混进频率轴会出现焦斑平移错位。第二ifftshift 的作用是把 fft2 输出的零频分量移回数组中心漏掉它会让传播相位 H 错位半个网格焦面强度图出现棋盘状混叠。第三kz 中的虚数项对应倏逝波应直接置零而不是写成 exp(−|kz|z)否则会引入虚假吸收。z 为负时同一段代码就是伴随传播。酉传播矩阵的伴随正是把 kz 换负号等价于反向传播。网格间距 dx 建议取 λ/(4~6)λ633 nm 对应 0.12~0.16 μmdx 太大则高频空间频率采样不足dx 太小则 FFT 内存按 N² 增长512×512 双精度复数阵已经需要几十 MB再往上要关注内存占用。2.3 单元相位响应的表格化FDTD/RCWA 跑一次优化里全用插值逆向设计最终要落到可制造的单元结构不能只算连续相位板。优化迭代时如果每一步都调用 FDTD/RCWA 计算单元响应几百个单元、几十个频率、几百次迭代算力上撑不住。常见做法是预先对单元库做一次扫描把柱径或边长、高度作为变量仿真透射相位和振幅存成表格优化循环里用插值代替仿真。% 单元库扫描结果示例来自FDTD, 数据量通常几十到几百条) d_list 80:20:400; % 纳米柱直径, nm phase_list [-1.2 -1.8 -2.5 -3.1 -3.8 -4.6 -5.2 ...]; % 透射相位, rad amp_list [0.93 0.92 0.90 0.88 0.85 0.82 0.80 ...]; % 振幅透过率 phase_lookup (d) interp1(d_list, phase_list, d, spline); amp_lookup (d) interp1(d_list, amp_list, d, spline);两个容易翻车的地方。一是相位数据必须先 unwrap扫描得到的相位通常卷绕在 [−π, π]直接插值会在卷绕点产生跳变放进梯度计算后直接崩掉保存表格前对相位做 unwrap插值时若柱径超出扫描范围返回 NaN 或截断不要外推。二是同一柱径在不同偏振、不同入射角下响应不同全偏振设计至少要准备 TE/TM 两组表插值时按通道分别查表不能共用一组数据。3. 伴随梯度引导的逆向优化MATLAB核心算法实现3.1 目标函数怎么定最大焦斑能量还是贴合目标强度分布目标函数决定最终光学效果的上限。常见两种写法。第一种是聚焦效率型直接最大化焦面中心区域光能量J −∫ mask(x)·|E_f(x)|² dx负号是因为优化器一般做最小化。掩模 mask 取以理想焦点为中心的小圆盘半径约 2~3 个艾里斑半径即 1.22λf/(2R) 的几倍。第二种是强度分布匹配型要求焦斑强度分布逼近指定曲线用于光束整形类器件J ∫ | |E_f|² − I_target |² dx两种写法结构不同梯度公式类似。聚焦透镜我一般用第一种再加一个漏光惩罚项避免能量被堆到焦斑外整形器件用第二种。mask zeros(Ny, Nx); [y, x] meshgrid((1:Ny)-Ny/2, (1:Nx)-Nx/2); mask(x.^2 y.^2 (focus_r/dx)^2) 1; % focus_r 约 2~3 个艾里斑半径mask 用硬掩模配合边界过渡带即可不需要高斯模糊。漏光惩罚项写成 β·∫(1−mask)·|E_f|² dxβ 取 0.05~0.2太大焦斑变松太小能量泄漏到外围。3.2 伴随梯度一次正向加一次反向拿到全部变量的梯度超表面透镜动辄几万到几十万个设计变量。对每个变量用有限差分求梯度一次迭代就是几万次前向传播伴随法只需要正向一次、反向一次梯度计算成本与变量数无关。推导路径不复杂。记透镜面出射场为 t_field amp·exp(jφ)焦面场 E_f P(t_field)P 是传播算子。损失 J −∫mask|E_f|²对 φ 求导用伴随算子 P† 把掩模误差反向传播回透镜面function [loss, dphi] compute_gradient(phi, amp, lambda, f, dx, mask) % 核心: 计算损失值和伴随梯度 t_field amp .* exp(1j*phi); % 透镜面出射场 E_f prop_asm(t_field, f, lambda, dx); % 前向传播到焦面 loss -sum(mask(:).*abs(E_f(:)).^2); % 焦斑能量损失 dL_dEf -2 * mask .* E_f; % 损失对焦面场的导数 E_back prop_asm(dL_dEf, -f, lambda, dx); % 伴随反向传播 dphi -2 * imag( conj(E_back) .* t_field ); % 相位梯度 end代码里 t_field 是每个像素的出射复振幅E_back 是焦面误差反向传播回透镜面的伴随场。dphi 的物理含义是每个像素的相位增加一个无穷小量时损失的变化率。这个式子的成立依赖传播算子的酉性角谱法恰好满足反向传播直接调 prop_asm 并把 z 取负即可。最容易错的是符号。这里损失定义为负焦斑能量所以 dphi 带负号如果读者把损失定义成正能量梯度符号必须反过来否则会不断推着相位走向全口径幅值振荡。调试时可以先用 64×64 小网格用有限差分核对解析梯度两边在小扰动下应一致到两位小数以上。3.3 主优化循环最快跑通的最小迭代框架有了梯度主循环本身很直接。先给梯度下降版本方便看清每一步在做什么。phi target_phase(r_grid, lambda, f); % 以理想相位为初始解 lr 0.5; % 相位更新步长, rad for iter 1:200 [loss, dphi] compute_gradient(phi, amp, lambda, f, dx, mask); dphi dphi / max(abs(dphi(:)), 1e-12); % 归一化, 防个别点发散 phi phi - lr * dphi; if mod(iter, 20) 0 fprintf(iter%d, loss%.4f\n, iter, loss); end end步长 lr 的单位是 rad。相位分布一次更新超过几个弧度容易绕进局部极值所以 0.1~0.5 rad 起步比较安全。梯度归一化是第二个关键超表面边缘单元和高频单元的梯度能差三到四个数量级不归一化收敛曲线会明显抖动。这个版本收敛偏慢实际项目我会换成 Adam 更新梯度公式不变仅把phi phi - lr*dphi换成动量更新。对超表面这种自由度极高、地形崎岖的问题Adam 头几十轮能快速压低损失后期再用学习率衰减微调。停止判据建议看焦面光强峰值和聚焦效率是否连续 20 轮变化小于 0.1%不要只看损失绝对值因为掩模尺寸直接影响损失的数值范围。循环里每轮记录焦点峰值和 η_focus这两个量更接近物理目标是否达成。3.4 从连续相位到纳米柱几何后处理映射连续相位优化结果不能直接加工要映射到单元库的离散柱径。最直接做法是逐像素找相位误差最小的柱径d_opt zeros(size(phi)); for k 1:numel(phi) phase_err abs(angle(exp(1j*(phase_lookup(d_list) - phi(k))))); [~, idx] min(phase_err); d_opt(k) d_list(idx); end这里用angle(exp(1j*(...)))把相位差包到 [−π, π]算的是循环相位差而不是线性差否则 ±π 交界处的柱径映射会出错。映射完成后把连续相位对应的焦场与离散柱径对应的焦场各做一次前向传播对比效率估算映射损失。若损失超过 5%多半是单元库相位密度不够或柱径扫描范围太窄先补单元库再重跑优化不要在装配端返工。4. 参数怎么调从“能跑”到“跑出好结果”4.1 漏光惩罚权重与效率指标漏光惩罚系数直接决定能量是被聚焦还是被摊开。系数太小优化器会把能量往外推系数太大焦斑扩展。我的参考区间是 0.05~0.15。先跑 20 轮看损失曲线如果焦面能量在下降但焦点邻域能量占比低于 30%把惩罚系数往下压如果能量分布太散不聚焦往上加。配一个效率统计函数来看状态eta_focus sum(mask(:).*abs(E_f(:)).^2) / sum(abs(E_f(:)).^2);η_focus 不是 Strehl 比但能快速提示能量是集中还是发散。经验上 η_focus 到 0.6 以上的设计才算及格。4.2 初始化影响从随机相位起步基本是浪费时间超表面逆向优化的损失地形非凸从随机相位初始化容易陷进局部极值损失不再下降焦点边缘却布满细小旁瓣。常见做法是用理想双曲相位做初值只加微小随机扰动打破旋转对称性phi target_phase(r_grid, lambda, f) 0.05*randn(Ny,Nx);0.05 rad 的扰动足够。对目标相位接近物理最优的透镜设计这比零相位起步收敛快几十轮。如果做的是像差校正、特定入射角或波长补偿先用几何光学推一版初始相位再进优化效果比纯随机好得多。4.3 学习率与 Adam 自适应更新梯度下降适合快速验证正式设计我一般用 Adam。MATLAB 里手写轻量版本并不复杂m zeros(size(phi)); v zeros(size(phi)); beta1 0.9; beta2 0.999; eps 1e-8; for iter 1:300 [loss, dphi] compute_gradient(phi, amp, lambda, f, dx, mask); m beta1*m (1-beta1)*dphi; v beta2*v (1-beta2)*dphi.^2; mhat m / (1-beta1^iter); vhat v / (1-beta2^iter); phi phi - lr * mhat ./ (sqrt(vhat) eps); end这里的 dphi 已经按全局最大值归一化过。学习率 lr 从 0.1 起步100 轮左右乘以 0.3200 轮左右再乘一次 0.3。Adam 的动量项会让相位在局部极值边缘冲出去对超表面不完全是坏事偶尔能越过小尖峰。4.4 制造可行性约束别让优化结果超出可加工范围纯相位连续空间里优化结果可能出现过窄柱间距或极端占空比实际工艺做不出来。我常用的低成本方案是两步交替先算纯相位再把结果投影到可行柱径集合最后用投影结果重新初始化继续优化。% 投影到可行柱径范围 d_mapped min(max(d_opt, d_min), d_max); phi_proj phase_lookup(d_mapped); % 用投影相位重新初始化下一轮优化 phi phi_proj 0.02*randn(Ny,Nx);核心思想是不要一步硬投影。硬投影一次就固定设计损失几乎必然回升把投影结果作为新初始解再优化让损失在制造约束边界上继续往下走最终收敛位置才同时满足物理和工艺要求。5. 逆向设计避坑指南五个必踩的坑5.1 现象优化后的相位边缘出现奇怪突变环原因口径边缘单元的前向传播贡献被硬截断角谱法在截断处产生衍射伪影梯度在边缘被异常放大优化器顺势在边缘堆积了无物理意义的大相位变化。 解决把口径掩模从硬边改成 2~3 个网格的渐变过渡带同时对梯度做边缘衰减。我最初做单波长设计时在边缘踩过一次光斑中心挺亮但旁瓣结构很不干净后来检查相位分布才看到这一圈“光环”。5.2 现象优化出来的焦点很亮效率数据完全对不上原因漏光惩罚权重太低或无惩罚能量堆进焦点邻域但在透镜口径内、焦点区域外分布着一圈环形能量η_focus 很低。 解决把 β 上调到 0.1 以上并且每轮打印 η_focus。看损失曲线会骗人看焦点峰值也会骗人只有 η_focus 这个比值能同时反映能量是否真的收拢。5.3 现象损失曲线像心电图跳上跳下无法收敛原因学习率过大或梯度没有归一化相位在某个像素上单步跨过 π。角谱传播对相位跳变很敏感一个像素的跳变就能让焦面场完全变样。 解决lr 降到 0.1 以下每次迭代前检查 max(abs(dphi))如果最大值超过几十 rad按全局最大值归一化而不是按像素自己做截断。遇到这类问题先看梯度量级再调迭代不要盲目加迭代轮数。5.4 现象换个随机种子优化结果完全不同原因非凸问题局部最优确实多尤其从随机相位起步时每个种子收敛到不同局部极值。 解决用理想相位初始化并行跑 3~5 个种子各加不同扰动选 η_focus 最高的结果。优化只跑一次不叫逆向设计。多种子并行成本不高对结果稳定性的帮助却非常明显。5.5 现象单元库映射后实测透镜和算法预期严重偏差原因单元响应表是孤立单元条件下仿真的实际阵列中相邻柱子的近场耦合会改变局部相位柱径差异越大偏移越明显。 解决在算法端控制设计时尽量避免相邻柱尺寸突变验证阶段用 FDTD 重建整个阵列的局部孔径做整阵仿真与孤立单元表结果对比。偏差超过容忍范围时把整阵仿真得到的相位响应送回去做一轮修正而不是硬用原结果加工。6. 验证方法从优化结果走向可交付设计6.1 设计完成后必做的三个验证第一步用映射后的离散结构复振幅不是连续相位做一次前向传播与连续相位结果对比 Strehl 比和焦斑形态。这一步完全不碰电磁仿真只是把“连续相位板”与“实际超表面”之间的差异量化出来成本几乎为零。第二步用 FDTD/RCWA 对一个小孔径阵列做整阵仿真截取中心区域验证单元耦合影响。如果焦点偏移超过一个衍射极限说明单元库数据有问题优化端要补相位补偿。第三步才是加工实测重点测 NA、聚焦效率和工作带宽。逆向设计最终好不好看制造后的透镜是否把大部分能量送到预定焦点而不是只看仿真数字。6.2 从单波长走向多波长损失函数合并的写法要扩展到多波长消色差把损失函数合并成多个波长项的加权和即可。每个波长独立前向传播、独立算伴随梯度梯度加权平均后更新相位loss 0; dphi 0; for wl lambda_list [l, d] compute_gradient(phi, amp_wl, wl, f, dx, mask); loss loss weight_wl * l; dphi dphi weight_wl * d; end注意每个波长的单元响应表要单独准备相位初始化也要考虑工作波段中心。这段代码把单波长框架几乎无缝扩展成宽带设计是我给团队做模板时的通用写法。6.3 我的最后建议接触过的优化迭代问题里超表面逆向设计算是最容易“主动加戏”的给损失函数加太多项、加过于复杂的约束、或者对初始相位做太多人工干预出来的仿真指标很漂亮加工后却经常翻车。我现在保持目标函数最少先跑通单波长聚焦设计确认“前向传播—伴随梯度—离散映射—整阵验证”这条链路闭环再加约束、加多波长、加视场。如果你正在搭这条链路我的建议是从 64×64 网格起步用理想相位初始化先把焦点和效率打印出来再放大到 256×256。迭代速度重要但能否稳定收敛更重要。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
