简介：《光学工程：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; z<0 表示反向传播（伴随用） % 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 μm；dx 太大则高频空间频率采样不足，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.3，200 轮左右再乘一次 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。迭代速度重要，但能否稳定收敛更重要。希望帮到你。

本文还有配套的精品资源，点击获取