超表面全息仿真全流程：从GS算法到CST/FDTD重建像

超表面全息成像这几年热度一直在线从显示、信息加密到AR光学随便一搜都是一堆论文。但真正上手做过仿真的人都知道论文里那张漂亮的“重建图像”背后是一整套跨工具链路先用MATLAB跑GS算法把目标图变成相位分布再把相位映射成超表面单元的几何参数最后丢进CST里用FDTD/时域求解器做全波仿真拿到接近真实的电磁响应。这三步每一步都有独立门道而最坑的恰恰是步骤之间的衔接——MATLAB算出来的相位怎么导进CSTCST里的边界、激励、监视器怎么设置才能对上焦我最初做这个方向时光是让GS算法收敛到一张干净的全息图就折腾了两周后来在CST里建阵列又踩了一堆边界条件的坑。这篇文章把我从“目标图像”到“仿真重建像”的完整做法、代码和踩坑记录都摊开讲适合正在做超表面全息仿真、相位恢复设计或者想复现论文结果的同学参考。1. 从目标图到重建像整条仿真链路先搭起来1.1 为什么说超表面全息的核心是相位工程传统的全息术记录的是物光的振幅和相位信息但普通感光材料只能记录强度所以经典全息要靠参考光干涉把相位信息编码成干涉条纹。计算全息把这个思路搬到了计算机里既然全息面的复振幅分布可以计算出来那就直接设计一个纯相位分布让它在特定距离处衍射重建出目标图像。超表面在这里扮演的角色就是一块能够精细控制空间各点相位的“人工界面”。它由亚波长单元构成每个单元的几何尺寸决定了对入射光的相位响应。你把这上千上万个单元的相位排布设计好入射光经过超表面后就会按设计好的波前传播在目标平面上相干叠加形成预设的图像。这里有个关键认知对于人眼观察的图像重建来说相位才是主角振幅扰动反而可以接受。GS算法要解决的核心问题就是“只知道目标图像的强度如何反推出全息面上每一点的相位”。1.2 从GS算法到CST/FDTD的四步链路整套仿真工作流可以拆成四步每步的输入输出非常清晰步骤工具输入输出关键点1. 生成全息图MATLAB目标灰度图像相位分布矩阵GS算法迭代获得纯相位全息图2. 映射单元参数MATLAB CST单元库相位分布矩阵超表面阵列几何参数量化相位查找表中匹配3. 全波仿真CST时域求解器/FDTD超表面三维模型近远场电磁数据边界条件、激励、监视器正确设置4. 提取重建像MATLAB后处理E场分布重建图像强度分布与理论重建对比很多初学者拿到GS算法就跑MATLAB里生成了漂亮的相位图然后直接卡在“怎么把相位图变成CST模型”上。实际上第二步的关键是建立一个“相位库”查表连续相位值量化成N个离散级别每个级别对应一个单元结构。第四步则是反过来CST算完电磁场后把监视器上的场强分布拉出来和MATLAB里的理论重建图对比。整条链路没有一步是可以省略的任何一步的参数设置不对最终图像都会一塌糊涂。2. GS算法一张相位图是怎么迭代出来的2.1 把远场衍射看成一次傅里叶变换GS算法Gerchberg-Saxton是1972年提出的相位恢复算法几十年了依然活跃在全息、光束整形、显微成像领域。它的核心物理基础很简单夫琅禾费衍射可以近似为傅里叶变换。想象一个平面光波通过全息面在很远的屏幕上形成衍射图案这个图案的复振幅分布就近似等于全息面复振幅的二维傅里叶变换。于是问题变成已知目标图像就是屏幕上的光强分布 求解全息面上的相位分布因为纯相位全息面振幅恒定这看起来很美但相位恢复是个病态问题没有闭式解。GS算法的解决办法是交替投影在全息面施加“振幅为1”的约束在目标面施加“振幅等于目标图像”的约束然后在这两个约束之间来回迭代。每一轮迭代都在修正相位振幅约束越来越接近目标最终相位分布收敛到一个能重建出目标图像的解。我常用一个生活化类比GS算法就像一个雕刻师手里有一块固定大小的木料纯相位全息面振幅恒定他不断把木料削成目标图案的形状目标振幅约束但每次削完看一眼反光形状全息面相位发现高光分布不对就再调整刀法最后反复打磨到形状和想要的大致吻合。2.2 交替投影迭代与收敛判据标准的GS迭代流程可以写成这样初始化相位矩阵通常用随机相位和目标振幅相乘构成初始场对该场做傅里叶变换得到全息面的复振幅分布保留全息面的相位把振幅强制设为1纯相位约束对修正后的全息面做逆傅里叶变换得到目标面的复振幅保留目标面的相位把振幅强制设为目标图像的振幅回到步骤2反复迭代直到误差不再下降关键问题是怎么判断收敛。我一般同时监控两个量归一化均方误差NMSE重建图像振幅与目标图像振幅之差这是最直接的指标。相关系数重建图像和目标图像的空间相关性相关系数高于0.9时视觉上就比较能接受了。实际操作中GS算法对初始相位很敏感。随机初始相位可能导致算法陷在局部极小重建图像有散斑噪声。我的习惯是同一目标图跑10次不同的随机种子取误差最小的一次。这10次计算也就几秒钟的事不要心疼这点算力最终效果的提升很值得。3. MATLAB实现GS算法代码、调参与误差统计3.1 可直接运行的GS脚本我给出一个我自己常用的GS算法脚本使用的是标准FFT快速傅里叶变换作为传播模型。代码注释尽量写清楚了每一行的作用你可以直接复制到MATLAB里跑。% GS.m - 标准GS算法生成纯相位全息图 % 用法: [phaseHolo, err] GS(target.png, 200, 0) function [phaseHolo, err] GS(imgPath, numIter, mode) % 读取目标图像并转灰度、归一化 tarImg imresize(rgb2gray(imread(imgPath)), [256, 256]); tarImg mat2gray(tarImg); % 目标振幅: 开方处理视觉上线性化 targetAmp sqrt(tarImg); N size(targetAmp, 1); % 初始场: 振幅用目标振幅相位用随机相位 rng(42); initPhase 2*pi*rand(N, N) - pi; field targetAmp .* exp(1i * initPhase); err zeros(numIter, 1); for k 1:numIter % 正向传播: 目标面 - 全息面 holo fftshift(fft2(fftshift(field))); holoPhase angle(holo); % 全息面约束: 纯相位振幅设为1 holoNew exp(1i * holoPhase); % 反向传播: 全息面 - 目标面 recon ifftshift(ifft2(ifftshift(holoNew))); reconAmp abs(recon); reconPhase angle(recon); % 计算重建误差 (考虑整体亮度比例因子alpha) alpha sum(sum(reconAmp .* targetAmp)) / sum(sum(targetAmp.^2)); err(k) sqrt(sum(sum((reconAmp - alpha * targetAmp).^2))) ... / sqrt(sum(sum((alpha * targetAmp).^2))); % 目标面约束: 保留相位替换振幅为目标振幅 field targetAmp .* exp(1i * reconPhase); % 前50次迭代每10次打印一次误差 if mod(k, 10) 0 fprintf(Iter %d, NMSE %.4f\n, k, err(k)); end end % 返回最后一次迭代的全息面相位 phaseHolo holoPhase; end跑完之后把相位矩阵phaseHolo保存下来后面步骤要用save(phaseHolo.mat, phaseHolo); % 顺便导出成灰度图查看 imwrite(mat2gray(phaseHolo), phaseHolo.png);这段脚本的核心不在于复杂而在于两个FFT之间的约束替换。这里有一个容易被忽略的细节fftshift和ifftshift一定要成对使用否则图像会整体平移错位。我最初写代码时漏了fftshift重建出来的图像中心图像跑到了四角排查了半天才发现是这种病。3.2 灰度图重建的两个关键技巧扩散器与加权GS如果你直接用上面的代码跑一张普通的灰度图比如人像第一次迭代之后你会发现重建图像对比度很低细节模糊。这是GS算法对灰度图天生的弱点原因在于灰度图在傅里叶频谱中能量主要集中在低频直接迭代难以把高频细节“挤”出来。解决这个问题有两个实用技巧。技巧一随机扩散器Random Diffuser在GS迭代开始前给目标振幅乘上一个随机相位模板diffuser exp(1i * 2*pi * rand(N, N)); field targetAmp .* diffuser;这个操作等于给目标图像人为加了一组随机相位让它在全息面的能量分布更均匀GS迭代的稳定性和收敛速度都会明显改善。重建后的光场会自动把这些随机相位“抹平”不影响最终肉眼观察到的强度图像。技巧二加权GSWeighted GS标准GS每次迭代都把重建振幅“硬”替换成目标振幅这样容易震荡。加权GS的思路是在替换振幅时做一个平滑混合w 0.7; % 权重一般0.5~0.8 field (targetAmp .^ w) .* (reconAmp .^ (1-w)) .* exp(1i * reconPhase);当w接近1时算法倾向快速逼近目标w偏小则迭代更“温和”不容易陷入局部极小。实际使用中我经常让权重随迭代次数动态变化前50次用w0.5之后逐步升到0.9这样兼顾收敛速度和重建质量。跑完GS后可以顺手画一下误差收敛曲线如果曲线下降后又有明显的反弹说明权重或者初始相位设置有问题换一个随机种子再试。4. 从相位图到CST超表面单元库设计与映射4.1 单元结构选型与相位覆盖MATLAB算出来的是连续相位的“理想值”现实中超表面单元只能提供有限的离散相位状态。所以在建CST模型之前必须想清楚用什么单元结构、怎么实现0到2π的相位覆盖。最常用的三类超表面单元介质纳米柱如非晶硅、氮化镓通过改变柱子的长和宽在透射模式下实现相位调控。优点是透射率高工艺成熟是目前可见光超表面的主流方案。等离子体V形天线/棒状天线通过局部表面等离激元共振调节相位结构简单但欧姆损耗较高更适合红外和太赫兹频段。几何相位单元Pancharatnam-Berry相位用各向异性单元如矩形柱对不同旋转角度实现2倍转角相位延迟。单元尺寸固定只需调节旋转角相位分布连续且相位跳变很小工程上很好用。我个人做可见光波段超表面全息最常选的是介质纳米长方柱加上几何相位方案。长方柱的长宽都可以作为变量能形成比较大的相位库再配合旋转角几乎是“怎么用都够”的灵活度。无论选哪种都要满足一个硬指标参数扫描范围内相位能覆盖完整的0到2π并且透射/反射振幅尽量高。如果相位覆盖不够重建图像上会有明显的周期性暗区或噪点。4.2 在CST中用参数扫描建立相位库假设选一个“硅纳米长方柱玻璃衬底”的结构周期p300nm柱高h600nm工作波长λ633nm。设计变量是长方柱的长Lx和宽Ly扫描范围取80nm到250nm。在CST里建相位库的步骤新建一个单元模型一块玻璃衬底加一个Si长方柱按实际尺寸画好。设置边界条件为unit cellX、Y方向周期边界Z方向开放。使用频域求解器Frequency Domain Solver设置波端口或平面波激励入射方向沿Z。把Lx和Ly设为参数运行参数扫描Parameter Sweep。提取固定频率点的S21透射系数的幅度和相位。扫描完成后整理成类似下面的查找表数据为示例编号Lx (nm)Ly (nm)透射相位 (deg)透射幅度111014000.922130160450.903150170950.8841701801400.8751902001850.8562102102300.8072302202800.7882502403300.75这个表就是GS相位映射的“字典”。注意不要只看相位幅度低的点要舍弃否则重建图像会有能量空洞。我在实际扫描中经常用保存相位和幅度到MATLAB脚本里的方式自动筛选只保留幅度高于某个阈值的参数组合例如0.7再从中挑选相位均匀分布的8个或16个点。4.3 量化映射与批量建模拿到相位库后把GS输出的连续相位量化成离散相位电平。量化方法很简单% phaseHolo: GS得到的连续相位范围 -pi ~ pi % lookupPhase: 从CST提取的离散相位向量 (弧度) [~, idx] min(abs(exp(1i*phaseHolo(:)) - exp(1i*lookupPhase.))); idx reshape(idx, size(phaseHolo)); % 每个idx对应一个单元几何参数组合我通常量化到8个或16个相位电平。8个电平的加工宽容度高、仿真量小重建图像质量也足够16个电平的图像噪声更低但单元种类多建模和加工都麻烦适合精细复现的场景。CST里批量建模的话可以用VBA宏或者通过CST的Brick组件脚本生成。核心思想是按相位索引查询几何参数然后在每个格点位置创建一个对应尺寸的长方柱。一个简单的VBA伪代码思路如下For i 0 To 19 20 x 20阵列 For j 0 To 19 idx PhaseMap(i, j) 从MATLAB导出的相位索引矩阵 Lx lookupLx(idx) Ly lookupLy(idx) x0 i * period y0 j * period Brick.Name pillar_ i _ j Brick.X0 x0 : Brick.X1 x0 Lx Brick.Y0 y0 : Brick.Y1 y0 Ly Brick.Z0 0 : Brick.Z1 h Next j Next i如果你不习惯VBA也可以从MATLAB里生成一个CSV表格列出每个单元的中心坐标和长宽然后在CST中按表格手动复制只是阵列一大效率就太低了。所以我还是推荐把VBA流程吃透超表面仿真建模的效率能高一个数量级。5. FDTD全波验证边界、激励、监视器的正确配置5.1 支撑CST时域求解器从单元库到有限阵列的身份切换CST的时域求解器基于有限积分技术FIT它和FDTD在离散方式和时间迭代思路上同源很多场景下大家直接把它称作FDTD仿真。用它来做超表面全息验证最大的优势是一次宽频计算就能拿到宽频带响应比频域求解器快不少。最容易翻车的地方是边界条件。做单元库提取相位时用的是unit cell周期边界模拟的是“无限大阵列”。但验证全息成像时目标图像是在一个有限大小的超表面阵列衍射下形成的如果继续用周期边界电磁波会无限重复重建图像完全不可见或者出现一堆重复的衍射级次。所以验证阶段要切换到开放边界X/Y/Z方向都设为Open (Add Space)或Open (Add Space PML)模拟自由空间。入射激励用平面波沿Z方向正入射电场极化方向沿X要和设计时的单元库激励一致。在超表面上方和下方预留足够空间加上PML吸收边界通常留出半波长到两个波长的距离。设置汇总可以参照下表参数设置说明求解器Time Domain SolverFIT/FDTD方法边界X/YOpen (Add Space PML)有限阵列不设周期边界ZOpen (Add Space PML)入射和透射方向开放激励Plane Wave正入射极化方向X网格最大网格步长约λ/12保证精度柱体内部至少2层网格采样频率单一频点或窄带取设计波长网格设置不要太保守太粗网格会让相位响应偏大重建图像模糊太细又会让仿真时间爆炸。以可见光633nm为例单元周期300nm用每波长12个网格的密度一个20×20阵列通常几十分钟能算完属于可以接受的范围。5.2 提取重建图像并对比理论结果全息成像的“像平面”在超表面上方的某个距离处。这个距离取决于你的成像架构傅里叶全息像平面在远场无穷远但实际可以用透镜把远场拉近到焦平面。菲涅尔全息像平面在近场距离超表面几十微米到几百微米CST仿真里可以直接放置监视器。对于CST仿真验证做菲涅尔型更直接给目标图像一个设计距离d例如100μm然后在这个位置添加一个E-field监视器场监视器Field Monitor监视频率设为设计波长对应的频率。仿真结束后导出该监视器上的电场强度分布在CST后处理中导出场的数据E-Monitor - Export - ASCII或MATLAB接口。用MATLAB读取矩阵取abs(Ex).^2 abs(Ey).^2总场强度就是重建图像的强度分布。把这个强度图归一化到0~1和前面GS算法理论重建出的强度分布并排对比。对比的时候不用追求像素级一致。CST是真正的三维全波电磁仿真包含了单元间耦合、损耗、网格离散误差等实际情况和理想FFT模型相比会有偏差。重点看三个指标图像轮廓是否可辨识如果重建图像还能看出目标图形的轮廓说明相位映射和仿真链路是通的。暗区是否出现异常的杂散光点如果出现大量亮点或周期性鬼像通常是相位量化不当或单元耦合导致的。整体亮度是否均匀亮度严重不均往往指向单元库中部分参数幅度太低。我第一次跑通时重建图像形状是对的但暗区有很刺眼的散斑噪声。查了半天发现是CST网格对柱子尖端部分剖分不足导致高频相位误差。把网格加密后噪声降了不少。这种问题不看仿真云图根本发现不了。6. 实测踩坑清单从仿真到复现论文的经验6.1 相位跳变是重建图像破损的头号元凶GS算法算出的相位分布通常是连续渐变的但经过量化后相邻相位值如果突然从接近0跳变到接近2π就会在超表面表面上形成一条“相位裂缝”。这种跳变直接对应单元结构参数的突变局部耦合和衍射效应会让图像上出现明显的条纹或暗线。缓解的方法有三个相位连续性约束在GS迭代中加入相位平滑约束或正则项使相位图在空间上变化更平缓。改用几何相位单元单元尺寸固定只旋转角度相邻单元的相位差就是旋转角差的2倍跳变量天然被控制。使用蓝噪声抖动量化量化时打乱误差分布避免均匀量化产生的周期性栅格图像背景噪声会更像“散斑”而不是“网格”观感好很多。6.2 单元间耦合让“查表法”失灵单元库是从“无限周期阵列”提取的也就是每个单元四周都是同样尺寸的结构。但当GS相位图里相邻单元的尺寸差异较大时实际阵列里每个单元的电磁环境都不是均匀周期环境单元间互耦会让透射相位偏离查表值。怎么发现这个问题我的做法是CST里抽一个2×2或者3×3的单元组分别用这个组内的相位组合去仿真对比单单元查表给出的透射相位的偏差。偏差超过30度就需要警惕了。应对手段有几种尽量让相位图在空间上平滑减少相邻单元的尺寸跳跃。选对耦合不敏感的单元类型比如高折射率介质柱在弱耦合条件下表现比金属等离激元单元好得多。如果条件允许把单元库从“尺寸维度”换成“旋转角度维度”的几何相位耦合一致性更好。6.3 仿真规模、网格与时间的三方权衡全息成像验证最常见的尴尬是GS算法生成的是512×512的相位图你不可能在CST里建几十万个单元做全波仿真。就算能建内存和求解时间也会直接让你怀疑人生。我的做法是“三步走”CST负责“取库”用单元仿真把离散相位库提取出来这一步是全波仿真精确度足够。MATLAB负责“合成”把提取的相位库数据当成每个单元的复透射系数用惠更斯原理或角谱法在MATLAB里做阵列级衍射计算得到理论重建图像。这个步骤几乎不花时间可以放心跑512甚至1024的阵列。CST负责“抽查”在整图中截取一个重点区域比如20×20单元用CST时域求解器做全波仿真对比这20×20区域上MATLAB合成结果的差异用来验证仿真的准确性。这套组合流程既保证了效率又保留了三维全波仿真的可信度。我在实际项目中几乎都是先用MATLAB快速遍历各种GS参数和单元库方案确定优化方向后再启动有限区域全波仿真省下大量算力。文章写到这里我还想多说一句超表面全息仿真的难点从来不是单个工具不会用而是跨工具的参数传递和物理概念转换。GS算法不是魔法它只是把“相位恢复问题”转化成“两个约束之间的迭代投影”CST时域求解器也不是黑箱把边界条件和监视器搞对了它就是一套标准的电磁场计算流程。只要把链路里每个环节的输入输出弄清楚你完全可以从零复现出一套可用的全息仿真结果。希望这篇记录能帮你少走一些我走过的弯路。