1. 为什么超透镜设计绕不开"脚本驱动全波仿真"

第一次把 CST 和 MATLAB 串起来跑超透镜参数扫描，我把 300 组单元胞丢进去就下班了。第二天回来，一半的 case 因为端口模式设置冲突直接报错退出，剩下一半的相位数据画出来是一条锯齿状的折线，根本没法用。那次之后我才明白，超透镜这类器件的仿真问题，本质不是"会不会用某个软件"，而是"能不能把几万次全波求解组织成一条可信的数据流水线"。CST 负责把麦克斯韦方程算准，MATLAB 负责把参数、几何、数据、判据全部管起来，两边各干各的活，中间靠一层薄薄的脚本粘合。这个组合看起来简单，真正让链路稳定可复现，要处理的细节远比想象中多。

这篇内容面向的是已经能用 CST 建单个纳米柱模型、也写过一点 MATLAB 脚本，但还没把整条链路跑通的人。我会按"单元胞相位库—目标相位布局—全口径验证—公差回代"的顺序拆开讲，重点放在每一步为什么这么做、哪里最容易翻车。文中参数以近红外通信波段 1550 nm 附近的硅基透射型超透镜为例，换成可见光波段的二氧化钛体系，思路完全平移，只是折射率和高度要重新算。

1.1 相位调控把设计问题变成了参数扫描问题

超透镜的物理内核很朴素：在口径上铺满亚波长尺度的纳米柱，每个柱子对入射光施加一个局域相位延迟，让整个口径出射的波前按设计好的形状弯曲。要让焦点落在距离 f 处，需要的相位分布是

φ(r) = -k(√(r² + f²) − f) + φ0，其中 k = 2π/λ

这个式子只回答了"每个位置需要多少相位"，没回答"什么几何能给出这个相位"。后者必须靠全波仿真去测。一个柱子的半径从 20 nm 扫到 200 nm，相位从 0 连续转到 2π 以上，中间还夹着透射率的起伏。这意味着你必须对几十到上百个几何组合逐个做电磁仿真，得到一张"几何—相位—透射率"的查找表。查找表的质量直接决定后面全口径仿真的成败，而这恰恰是纯手工操作最容易出错、最需要脚本化的环节。

更麻烦的是，超透镜不是单一结构的重复，口径上每个位置的柱子参数都不同。500 μm 口径、300 nm 周期，就是大约 170 万个单元。全部真实建模做全波仿真，今天的硬件也扛不住。所以实际流程一定是分层：单元胞用周期性边界条件算准单柱响应，全口径用某种降规模手段验证整体波前。这两层的网格、边界、激励方式完全不同，用错地方结果会离谱得很自然，你还看不出哪里错了。

1.2 CST 与 MATLAB 的分工边界在哪里

我见过两种极端做法。一种是全在 CST 里用 VBA 宏跑，参数存在 Excel 里读进来，结果也写回 Excel，能跑但改一个参数要动宏代码，调试体验很差。另一种是全在 MATLAB 里用自己写的 RCWA 或 FDTD，灵活但精度和材料色散模型要自己搭，遇到复杂结构就力不从心。比较务实的做法是让 CST 只做它最擅长的事：给定几何和材料，把 S 参数或场分布算准；其余的参数生成、批量调度、数据清洗、插值查表、后处理成像全部放 MATLAB。

这样切分的好处是，仿真内核和设计逻辑解耦。换材料、换波长、换口径形状，只需要改 MATLAB 侧的目标相位函数，CST 模型文件保持不变，通过参数接口传值就行。反过来，如果哪天要换求解器（比如单元胞改用别的工具算），MATLAB 侧只要保证输出同样格式的相位库，后面全流程不受影响。

接口层面，Windows 下最直接的是 COM 自动化。CST Studio Suite 注册了CSTStudio.Application这个 ProgID，MATLAB 用actxserver拉起来，拿到主对象后就能打开工程、写参数、触发重建、启动求解器、读结果树。这条路的好处是不用装额外的工具箱，坏处是 COM 调用有开销，一次参数写一次重建，几千次循环里如果每次都Rebuild整个模型，光重建就能耗掉几小时。所以脚本设计上要做批处理：把能合并的参数变化合并，能用一次重建解决的不要拆成多次。

另一条路是 CST 自带的 MATLAB 接口（新版里叫 CST Studio Suite 的 MATLAB 集成），封装程度更高，语法更接近面向对象，但版本之间方法名有差异，跨版本迁移时要小心。我的建议是先用 COM 把流程跑通，摸清每一步到底发生了什么，再去用封装更高的接口。

2. 把单元胞仿真做成可批量调用的"相位库工厂"

单元胞模型是整个流程的地基。它的任务只有一个：给定一组几何参数，返回该结构对正入射平面波的透射复振幅。听起来简单，但模型里任何一个设置出错，返回的相位都会系统性偏移，而且这个偏移在后面的全口径仿真里会被原样继承，很难追溯。

2.1 单元胞模型里最容易被忽略的三处设置

第一处是边界条件。单元胞模拟的是无限周期阵列中的一个柱子，所以 X、Y 方向必须设成 unit cell（周期性）边界，Z 方向用 Floquet 端口。很多人图省事把 X、Y 设成 open 加对称面，结果相当于在仿真一个孤立柱子，相邻柱子的近场耦合完全丢失，相位会偏几十度。在周期 p 接近 λ/2 的时候，这个偏差足以让聚焦效率掉一半。

第二处是端口到结构的距离。Floquet 端口不能贴着柱子放，否则倏逝波没衰减完就被端口吸收，S 参数会带上伪影。经验值是端口面距离结构上下表面各留 λ/2 以上，如果结构高度较大或者有高折射率基底，留到 0.7λ 更稳妥。这个距离会影响 S21 的绝对相位，但只要每次仿真保持一致，相对相位分布不受影响，所以关键是"固定"，不是"最短"。

第三处是材料色散模型。硅在 1550 nm 附近折射率约 3.48，但虚部不为零，有微弱吸收；如果用常数折射率模型，透射率会被高估，相位倒是影响不大。二氧化钛在可见光波段相对透明，但不同制备工艺的薄膜折射率差异能达到 0.1 以上，直接影响所需柱高。我的做法是把材料数据单独存一份文本文件，MATLAB 读进来后通过参数接口写进 CST 的材料定义，而不是在 CST 界面里手填。这样换材料、换数据源时只改一处。

2.2 Floquet 端口与相位解缠的实操细节

亚波长周期下，只有 0 阶衍射是传播模，其余都是倏逝模。所以端口模式数设成两个（TE 和 TM 的 0 阶）就够了，多设只会拖慢求解。正入射时 TE 和 TM 简并，取其中一个的相位即可；如果要考虑斜入射或者偏振转换，两个都要。

单频点求解用频域求解器更划算。时域求解器一次能出宽带结果，但你只需要一个波长，频域求解器在这个场景下的速度快好几倍，而且 Floquet 端口的相位提取更干净。这一点在我早期的流程里被忽略了很久，白白多花了不少机时。

读出来的 S21 相位是包裹在 (−π, π] 区间里的，必须解缠才能看到真实的相位累积。MATLAB 里一行unwrap(angle(S21))看似解决，但要注意扫描顺序：如果半径是从小到大扫的，相位随半径单调变化，unwrap没问题；如果扫描顺序乱了，或者参数空间是二维的，就必须先排序再解缠。我吃过这个亏——一次二维扫描（半径 × 高度）忘了排序，unwrap在数据跳变处凭空加了 2π，导致整张查找表错位。

% 先把扫描结果按半径排序，再解缠 [~, idx] = sort(scanTable.R); scanTable = scanTable(idx, :); scanTable.phase = unwrap(angle(scanTable.S21));

还有个细节：相位参考面。S21 的相位包含了波在端口之间传播的路径相位。如果你在扫描过程中改变了柱高，参考面到结构表面的距离会变，路径相位也跟着变，看起来像是柱子本身的相位在变，其实混入了几何路径项。解决办法是在后处理里减掉k0 * (端口间距 − H)这一项，或者干脆固定端口位置、只改结构高度，让路径项成为可计算的常数。

2.3 MATLAB 驱动脚本的骨架

下面这段是我现在用的批量扫描骨架，思路是"一次连接、批量改参、批量求解、批量读结果"，而不是每次循环都开关 CST。COM 对象建立和销毁的开销在几百次循环里能占到总时间的 20% 以上。

cstApp = actxserver('CSTStudio.Application'); mws = cstApp.OpenFile(fullfile(workDir, 'unit_cell.cst')); resTree = invoke(mws, 'Resulttree'); nCase = numel(scanTable.R); for k = 1:nCase invoke(mws, 'StoreDoubleParameter', 'R', scanTable.R(k)); invoke(mws, 'StoreDoubleParameter', 'H', scanTable.H(k)); invoke(mws, 'Rebuild'); solver = invoke(mws, 'Solver'); invoke(solver, 'Start'); obj = invoke(resTree, 'GetResultFromTreeItem', ... '1D Results\S-Parameters\SZmax(1),Zmax(1)', 'oh'); x = invoke(obj, 'GetArray', 'x'); y = invoke(obj, 'GetArray', 'y'); scanTable.S21(k) = y(1) + 1i * y(2); % 具体格式看版本，先看一眼返回 end

方法名在不同版本里有出入，GetResultFromTreeItem的第三个参数（结果分量名）尤其容易写错。跑之前先在 MATLAB 里methodsview(mws)扫一眼可用方法，比翻文档快。返回对象如果是复数数组，注意它是分成实部虚部两个数组返回的，还是逐元素交错，这个必须实测确认，否则相位会整体错乱。

3. 从相位库到全口径布局：查找表的建立与陷阱

相位库建立起来后，剩下的事看起来是纯粹的插值查表：给定每个位置的目标相位，找到最接近的几何参数，然后把上百万个柱子铺到口径上。真正做过的人都知道，这里的问题不在于算法复杂度，而在于离散化带来的取舍。

3.1 目标相位面的几种常见形式

聚焦是最基本的形式，前面给过公式。实现时要注意把相位对 2π 取模后再查表，因为相位库本身只覆盖一个 2π 周期（有的结构能覆盖 2π 以上，多出来的部分其实是重复的）。取模之后目标相位落在 [0, 2π) 区间，和库的区间对齐。

偏折的形式是 φ(x) = −k x sinθ，只在一个方向上线性变化。涡旋光束是 φ(θ) = lθ，l 是拓扑荷数，相位绕中心一圈累积 2πl。这三种可以叠加，做出聚焦加偏折、聚焦加涡旋的组合功能。叠加时直接相加再取模，物理上对应波前的矢量叠加，是对的。

有一条容易被忽视：如果基底不是空气，焦点位置要按基底折射率修正。光从基底侧出射，焦距在介质里，公式里的 f 要乘上 n_sub，否则焦点位置会系统性偏近。我第一次做基底侧聚焦的时候没注意，仿真出来的焦点位置和设计差了将近 30%，排查了半天才想起来这一步。

3.2 相位跳变与透射率的取舍

扫描半径时，相位不是严格单调的。典型情况是：半径增大，相位先缓慢增长，到某个点附近因为模式共振相位快速跳变，然后继续增长，中间可能出现局部回环。这种非单调性会让同一个目标相位对应多个半径，选哪个？

判据是透射率。回环区域往往伴随透射率下降（共振吸收或强烈反射），所以选透射率高的那一支。具体做法是把查找表按相位排序，对重复的相位区间按透射率降序保留，再做单调化处理。这样得到的库虽然损失了一点相位分辨率，但整体效率更高。

[tbl, idx] = sortrows(scanTable, 'phase'); [tbl, ia] = unique(round(tbl.phase * 1e4) / 1e4, 'first'); % 相位重复处取首个 tbl = tbl(ia, :);

另一个取舍是相位量化。理论上每个位置都应该用精确匹配的半径，实际操作中会把相位量化成 8 档、16 档或 32 档，同一档的柱子完全相同。量化档数越多，波前越接近理想，但建模对象数和仿真时间线性上升。实测下来，8 档量化的聚焦效率比连续分布低 5% 到 8%，16 档能压到 2% 以内，32 档再往上收益很小。所以我的默认配置是 16 档，做最终验证时用 32 档对比一下确认收敛。

量化还有一个副作用：焦斑旁瓣会升高。因为量化相当于在相位面上加了周期性误差，误差的空间频率和量化台阶数相关，会在远场产生额外的衍射级。如果应用对旁瓣敏感（比如成像），档数要往高了取。

3.3 PB 相位与传播相位的混用边界

传播相位靠改变柱子的横向尺寸来改变等效折射率，从而改变相位累积。它的优点是 broadband，缺点是尺寸变化范围大，相邻柱子之间容易发生耦合，而且高度固定时相位覆盖范围有限。

PB 相位（几何相位）靠旋转柱子的朝向实现，相位等于旋转角的两倍（圆偏振入射时），与尺寸无关。优点是相位覆盖天然是 2π，色散小；缺点是必须用圆偏振光，而且偏振转换效率不是 100%，一部分能量会跑到正交偏振上。

两者可以混用：用传播相位提供基础相位，用 PB 相位做精细调节。混用的前提是单元胞的偏振转换特性要单独扫描一遍，把转换效率作为第二个筛选维度。如果转换效率低于 70%，混用的收益会被损耗吃掉。我一般只在需要宽带工作或者大数值孔径的场景才考虑混用，窄带单波长应用纯传播相位更省事。

4. 全结构验证仿真：网格、边界与效率口径

单元胞算得再准，也只是局部响应。全口径仿真要回答的是：这些柱子拼在一起，焦点的位置对不对，焦斑大小符不符合衍射极限，效率到底有多少。这一步的规模比单元胞大两三个数量级，设置上的任何疏忽都会被放大。

4.1 从单元胞到全口径，网格策略完全不同

单元胞用周期性边界，网格只要在柱子附近加密就行，全局可以用较粗的六面体网格。全口径结构用开放边界加平面波激励，网格必须覆盖整个口径和焦区，总量直接爆掉。

我的做法是三层网格：全局网格按 λ/12 设定，保证传播波的基本精度；在纳米柱所在的薄层加局部网格，密度提到 λ/30 左右，因为高折射率柱子的场变化剧烈；在焦平面附近再加一层中等密度的网格，保证焦点区域的场采样足够细。三层网格用 CST 的 local mesh 功能分别指定，不要试图用一套全局网格通吃。

网格收敛性必须验证。做法是把全局网格加密一档（比如从 λ/12 到 λ/18），看焦点位置和焦斑宽度变化多少。如果变化小于 3%，认为收敛；如果超了，说明网格还不够，或者某处有几何奇异点没处理好。柱子的侧壁如果是完美的垂直圆柱，网格在圆形边界上的阶梯近似会引入误差，加密能缓解，换更细的局部网格更有效。

4.2 焦平面后处理的三个口径问题

仿真跑完，拿到焦平面的场分布，接下来的计算有三个地方容易出定义分歧。

第一个是焦点位置的判定。取场强最大值点？还是取光轴上强度最大值？两者在大数值孔径下会不一样，因为焦斑可能不是完美的圆形对称。做轴上聚焦时用光轴最大值更符合物理定义；做离轴聚焦或者有像差时，用全场最大值更合理。报告结果时要把用的哪个说清楚，否则和别人对比会打架。

第二个是焦斑宽度的定义。半高全宽（FWHM）是最常用的，但对艾里斑形状的焦斑，主瓣两侧很快降到零然后回升到旁瓣，半高位置可能落在主瓣边缘之外，导致 FWHM 偏大。稳妥的做法是只在一阶暗环以内找半高位置，超过就认为定义失效。

第三个是效率的分母。透射效率的分母是入射到口径上的总功率，聚焦效率的分母也是它，分子是焦斑主瓣内的功率。但"主瓣内"的积分窗口取多大，结果能差 10% 以上。我一般取到一阶暗环位置，同时在报告里注明窗口大小，方便复现。

% 焦平面强度积分，窗口取一阶暗环 I = abs(E).^2; mask = rho <= rFirstDark; P_focus = sum(I(mask)) * dA; eta = P_focus / (I0 * A_aperture);

4.3 内存不够时的降规模手段

500 μm 口径在 1550 nm 下是 300 多个波长，全波直接算基本不现实。硬件允许的话，把口径缩到 100 μm 先验证波前形状，趋势对了再考虑放大。缩口径会改变数值孔径和焦斑，但相位分布的正确性能在缩比模型上看出来。

第二个手段是对称性。平面波正入射、结构轴对称时，可以用两个对称面把计算域压到四分之一。CST 里设成 electric 和 magnetic 对称面组合，具体哪个方向用哪种，取决于入射场的偏振方向，设错了会把模式滤掉。这个技巧在单元胞上不明显，在全口径上能省 75% 的内存。

第三个手段是相位量化加结构复用。16 档量化后，口径上只有 16 种不同的柱子。虽然它们的空间位置不同，不能简单阵列复制，但可以把每种柱子做成一个参数化组件，用脚本批量放置，同时共享网格模板。有些版本的 CST 对大量重复对象的网格生成有优化，对象数从百万降到万级以后，前处理时间能从几小时压到几分钟。

5. 踩过的坑：从脚本报错到结果不收敛的排查链路

下面这几个问题我都实际遇到过，排查过程比结论更有参考价值。

5.1 COM 连接失败的排查顺序

症状是actxserver报错，或者连接成功但OpenFile失败。排查按这个顺序走：先确认 CST 的位数和 MATLAB 的位数一致，64 位 MATLAB 连 32 位 CST 必然失败，这个坑在装了多个版本的工作站上特别常见。再看 CST 是否已经以交互模式开着，有时候手动打开的实例会占用 COM 接口，脚本再连会拿到一个空句柄。第三看路径里有没有中文或空格，COM 传路径对特殊字符的处理不稳，最稳妥是先把工程拷到纯英文短路径下。

连接成功但求解器不启动，多半是 license 被占满或者求解器类型没指定。脚本里Solver.Start之前要确认求解器已经切换到频域，默认可能是时域，而这个切换在 COM 里不是靠改一个属性，而是要调用具体的方法，版本之间命名不同。

5.2 相位曲线出现锯齿的根因

前面提过，第一次跑出来的相位是锯齿状。拆开看有三个叠加原因。一是没解缠，包裹的相位直接画出来当然是锯齿。二是扫描顺序乱，unwrap在跳变处误判。三是端口距离太近，倏逝波污染让 S21 在某些频点出现虚假的相位抖动。

第三个原因最隐蔽，因为它不表现为整体偏移，而是局部毛刺。判断方法是把端口距离从 0.3λ 加到 0.7λ 重跑一遍，如果毛刺消失，就是它。这个验证花不了多少时间，但能省掉后面几天的怀疑。

还有一种锯齿来自网格：柱子的圆形截面在六面体网格里是阶梯近似的，半径连续变化时，阶梯图案会突变，导致相位出现台阶状抖动。加密局部网格能压平，或者干脆把柱子建成方柱（精度略降但网格规整），看你更在意哪个。

5.3 效率虚高的常见误判

有一次算出聚焦效率 92%，高兴了半天，后来发现是分母算错了。入射功率密度是按平面波的 |E|²/2η 算的，但 CST 的平面波源有归一化设置，导出场的幅值单位可能不是 V/m 而是归一化值，直接套公式会差一个常数因子。

另一个误判来源是没有扣除基底反射。透射型超透镜的基底和空气界面本身有 4% 左右的菲涅耳反射，如果效率定义里包含这部分，永远到不了 100%。合理的做法是把基底另一侧的透射功率作为分母，只算超透镜结构本身的贡献。

第三个误判是焦斑窗口取得太大，把旁瓣甚至背景散射都算进主瓣。用一阶暗环做窗口边界，同时在报告里给出焦斑的完整剖面图，别人一眼就能看出窗口取得合不合理。

6. 公差与工艺回代：把仿真结果拉回现实

仿真跑出漂亮结果只是第一步。真正做出来之前，必须回答一个问题：工艺误差会让性能掉多少。这个问题不解决，设计再优也只是纸面数字。

6.1 参数化公差的批量跑法与快速估算

工艺上最常见的误差是柱子半径的偏差，光刻和刻蚀的综合误差在 ±5 nm 到 ±15 nm 之间。直接做法是蒙特卡洛：对每个柱子半径加随机扰动，重新生成全口径模型跑仿真。一次全波仿真几小时，跑几十次不现实。

更聪明的做法是利用相位库的局部斜率做快速估算。相位对半径的导数 dφ/dR 在库里已经有了（相邻扫描点的差分），半径偏差 ΔR 对应的相位误差就是 dφ/dR × ΔR。把这个相位误差叠加到目标相位面上，再用标量衍射积分算远场，几分钟就能得到统计结果。这个方法在相位误差不太大（小于 0.3 rad）时和全波结果吻合得很好，可以用来快速筛选工艺窗口。

dphi = gradient(tbl.phase, tbl.R); phiErr = dphi(idx) .* (randn(size(idx)) * sigmaR); phiTotal = phiTarget + phiErr; E_far = scalarDiffraction(phiTotal, lambda, z);

确定工艺窗口以后，再挑两三个极端 case 用全波验证一遍，确认快速估算没失效。这个"快速筛 + 少数验证"的组合，比全量蒙特卡洛省两个数量级的时间。

6.2 侧壁角与刻蚀深度的等效建模

实际刻蚀出来的柱子不是完美圆柱，侧壁有 85° 到 88° 的倾角，底部比顶部宽。严格建模要建圆台，网格量和对象复杂度都上升。等效做法是把圆台等效成直径取中值的圆柱，前提是倾角误差小于 3°，此时等效误差对相位的影响在 0.1 rad 以内，可以接受。倾角更大时，老老实实建圆台，或者用等效折射率法把锥形部分折算成一层渐变介质。

刻蚀深度误差比侧壁角影响更大。传播相位正比于柱高，5% 的高度偏差在 1550 nm、柱高 700 nm 的体系里对应约 0.3 rad 的相位误差，接近量化误差的量级。这个误差是全口径一致的（同一片晶圆上刻蚀深度基本均匀），所以它不会让波前畸变，只会让整体相位偏移，表现为焦点位置沿轴移动。补偿办法是在设计阶段把目标焦距按实测深度微调，或者在版图里加一层可调的余量。

这套流程我从单元胞模型搭起，前后改过三四版，最大的体会是：脚本的稳定性比单次仿真的精度更重要。一个能连续跑 500 个 case 不出错、结果格式统一的流水线，价值远高于一个手工调出来的完美单点。真正卡住项目进度的，从来不是求解器算得不够准，而是数据对不上、参数传错了、结果文件命名乱了这类琐碎问题。所以我现在写驱动脚本，第一件事是把日志和断点续跑做进去，每个 case 跑完立刻落盘，中途崩了能从最后一个成功的点接着来。这个习惯帮我省下的时间，比任何网格优化技巧都多。