数字型太赫兹带通滤波器的逆向设计及优化*

居学尉 张林烽 黄峰† 朱国锋 李淑锦 陈燕青 王嘉勋 钟舜聪 陈盈4)‡ 王向峰††

1) (福州大学机械工程及自动化学院,福建省太赫兹功能器件与智能传感重点实验室,福州 350108)

2) (CAD/CAM 福建省高校工程研究中心(莆田学院),莆田 351100)

3) (智能配电网装备福建省高校工程研究中心,福州 350108)

4) (福州大学物理与信息工程学院,福州 350108)

针对高性能太赫兹功能器件的规范化设计需求,本文将智能逆向设计方法应用于太赫兹带通滤波器的设计与优化中.建立与数字空间映射的亚波长超表面等效模型,从设定器件的目标功能和约束条件出发,利用智能算法探索整个解空间中的全部可能结构,迭代寻优至最优结构图案.本文利用搭建的逆向设计框架设计了中心频率为0.51 THz、带宽为41.5 GHz、插入损耗为–0.1071 dB 的太赫兹带通滤波器.与传统的人工正向设计相比,逆向设计方法可解构出窄带、低插入损耗、带外抑制强、极化稳定性强的带通滤波器.

1 引言

太赫兹(terahertz,THz)带通滤波器(bandpass filters,BPFs)是指允许通过某一特定频段的THz 波、同时将其他频段分量衰减到极低水平的滤波元件,因此广泛应用于THz 通讯、成像及隐身等领域[1–3].THz BPFs 的实现方式有很多种,包括全封闭金属结构波导[4]、光子晶体[5]、频率选择表面[6,7]、超材料[8]等.传统的设计方法依赖于设计者的经验和光学理论,在不改变主体结构的前提下,通过人为不断优化单元结构的特征参数以实现特定性能,设计流程复杂缺乏规范化,显然已不适合于未来光子器件智能化、精准化、设计周期短的发展趋势.

高效率、系统化的设计需求牵引着设计形式的更新与变革,智能逆向设计方法为高性能、多功能化的太赫兹器件设计提供了新的解决办法.该方法利用人工智能算法加速光学器件的设计和优化,不同于从器件结构向目标性能的正向设计方法,逆向设计方法是从器件的目标功能出发,通过人工智能算法探索整个设计区域的全部可能结构,自主迭代优化来寻找符合功能需求的结构.通过逆向设计方法设计出的光学器件具有非规则结构,但也能实现目标功能,可激发设计者的创新潜力并拓宽其设计思路.现如今人工智能算法广泛应用于光学器件的逆向设计中[9],大致可分为基于群体的启发式算法和基于训练集的机器学习类算法两种设计模式.如利用遗传算法反复迭代设计出了超材料吸波体[10–12]和纳米光子波长路由器[13]等,利用深度学习方法准确预测了三维手性超材料[14]和超表面天线[15]等,还有结合机器学习和进化算法逆向设计出了石墨烯超材料[16]和壳基机械超材料[17]等光学器件.

然而在THz 波段,逆向设计方法的研究还相对匮乏.本文基于两种智能算法开展了“数字型”(像素型) THz BPFs 结构的逆向求解和结构优化研究,智能算法依据设定的目标功能和约束条件在整个解空间中进行迭代寻优至最优结构图案.该方法不需要任何先知的结构或形状,仅需要指定设计区域及目标功能函数等参数,即可实现从性能设定到目标结构的逆向设计,智能逆向设计方法作为创新设计的有益补充不限于此案例,可促进现代光学器件设计从被动优化到逆向创造.

2 逆向设计方法

若要操纵太赫兹波完成更为复杂的功能,则需要在设计太赫兹功能器件时考虑更多的因素,而传统的正向设计方法已经无法满足其需求.逆向设计方法打破了传统正向设计方法多自由度受限的困境,以器件目标性能为导向,利用智能算法探索设计区域所有可能的结构,直到设计出符合目标性能的结构分布.该方法不仅节约人工成本和缩短设计周期,还让设计过程更加规范化和精准化.本文以光谱分析系统/成像中重要构成器件太赫兹带通滤波器为研究案例,开展数字型太赫兹功能器件的逆向设计研究,具体数字编码化仿真建模步骤如下.

2.1 基于直接二进制搜索算法的数字化仿真建模

直接二进制搜索(direct binary search,DBS)算法本质上是一种遍历搜索算法,通过迭代搜索整个设计区域从而得到最优全息图.假设设计区域为正方形,并按照一定的策略进行离散编码化,划分成n×n(n取决于设计精度)个方块(像素网格),该算法从第1 个像素点逐一更改图像中的像素点的逻辑状态(从0 变为1,或从1 变为0),在每次更改图像的一个像素点后都需要重新计算目标函数,即每次迭代都要进行n×n次运算.DBS 算法已经被广泛应用于偏振分束器和功率分配器的设计和优化中[18,19],当然该思想也可应用于THz BPFs 的逆向解构中,利用其对传统频率选择表面(FSS)的单元结构、特征尺寸和像素周期等参数进行优化,其目的是寻找性能最佳的介电常数分布,以实现设计性能指标.基于DBS 算法的数字化逆向仿真建模步骤如下.

1) 物理模型抽象为数学语言.传统的十字型FSS BPFs 的结构示意如图1(a)所示,FSS 的单元结构是二维周期结构中最小的元素(红虚线方框标注).而数字编码化是指对该最小元素继续进行切割,划分成为n×n个可寻址的像素点(红实线方框标注),其结构示意如图1(b)所示,与之相对应地产生一个长度为n2的编码串.每个像素点有两个逻辑状态,“0”代表空气(金属材料被刻蚀),“1”代表金属(金属材料没有被刻蚀),因此与物理模型一一对应的编码串成为二进制字符串,至此完成了亚波长结构与离散编码序列之间的映射建模.

图1 (a)传统的十字型FSS;(b)数字型THz BPFs 中像素点分布和1/8 对称结构示意图(红三角形虚线框所示)Fig.1.(a) A traditional cross-type FSS;(b) schematic diagram of the pixel distribution and 1/8 symmetrical unit structure surrounded by red triangle in a digital THz BPFs.

2) 设定器件理想性能的目标函数并定义品质因子(figure of merit,FOM).本文在设计THz FSS 时,目标是中心频率处透射率越高,带宽越窄,带外透射率幅值越小,因此该目标性能指标可近似定义为一维高斯函数,如图2(a)所示.在Matlab软件中可通过调整函数gaussmf(x,[AB])中的A/B参数控制其形状,参数A控制通带的半峰全宽(full width at half maxima,FWHM),参数B控制中心频率fgoal的位置.FOM 值是将结构的实际光学透射率曲线与目标透射率曲线进行比较,定义为多点频率处透射率与理想值差值的平方和,故本文FOM 值定义如下:

图2 (a)理想性能的目标函数;(b)逆向设计方法流程图Fig.2.(a) Ideal performance indicators;(b) the flow chart of reverse design method.

其中,T1为fgoal处的透射率;T2,T3分别为fgoal±BFWHM/2 处的透射率(BFWHM为设定的带宽);T4,T5,T6,T7为均匀分布在通带外各频率点处的透射率,如图2(a)所示.

3) 初始化并迭代寻优.初始化过程是通过在Matlab 数值编译软件中随机生成一串长度为n2并由0 和1 数字组成的二进制编码串,经解译后将等效的物理模型导入电磁场仿真软件(FDTD Solutions)中进行结构的光学响应(透射率)计算,返回的仿真结果代入(1)式中计算初始FOM 值.然后,从第1 个像素点逐一更改初始编码串中像素点的逻辑状态,令其状态0→1 或1→0,解译新的编码串并对新结构进行仿真计算,获得其FOM值.与先前的FOM 值进行比较,如果FOM 值得到改善,则保留更改编码串点处的逻辑状态,反之则保留原来的逻辑状态.在DBS 方法中,逐点轮询各像素单元的逻辑状态,总共需要进行n2次更改,称作一次迭代,并保留本次迭代中最优的FOM值.重复上述Matlab 与FDTD 联合优化步骤,经过多次迭代,当FOM 值收敛,性能不再提升时终止,最终得到满足设计要求的最优结构图案,迭代寻优流程图如图2(b)所示.

值得权衡的是像素点越小,输出的器件功能与设计性能越匹配,不过随之带来的是对应编码长度呈指数型增长,消耗巨大的计算资源.因此,采用旋转中心对称结构并引入旋转对称机制进行编码,使编码串长度减小为原有长度的1/8,其1/8 对称结构如图1(b)所示.该处理方式不仅大大减少了仿真计算时间,即每次迭代只需完成n(n+2)/8 次可寻址像素点的翻转和仿真计算;还使器件具有极化不敏感特性,可以直接抵消电场极化带来的影响,即器件结构不但关于x轴和y轴对称,还关于两条对角线对称[20].

本文器 件设计 参数如 下:fgoal=0.5 THz,BFWHM=0.05 THz,单元结构周期L=150/fgoal=300 μm,单元结构划分为20 × 20 个像素点,每个像素点边长l=L/20=15 μm,材料为铝箔,其等离子体频率为2.24×1016rad/s,阻尼系数为1.12×1014rad/s[21],厚度d=10 μm.FDTD 中其他参数设置如下: 仿真时长t=200 ps,平面波光源频率范围为0.1—1.5 THz,全局网格划分精度为5 μm,局部网格(作用区域)划分精度为1 μm.

2.2 基于离散粒子群优化算法的数字化仿真建模

运用DBS 算法优化设计区域时,会逐点轮询所有像素单元的逻辑状态,耗时较长且非常容易陷入局部最优解,为此采用全局搜索能力强、收敛速度快且精度高的离散粒子群优化(binary particle swarm optimization,BPSO)算法进行逆向设计[22].BPSO 算法是一种基于群体智能的进化算法,能根据环境中获得的适应度信息,应用某些类型的算子来更新个体,从而使种群中的个体向更优的解区域移动.因此,BPSO 算法善于在离散空间中寻找最优解,每个粒子代表可行解集中的一个解,定义某一时刻粒子i的位置X和速度V如下:

位置向量X(i)中每个分量元素xn(i)为0 或1,对应映射模型的数字编码串;速度向量V(i)中分量元素vn(i)为特定范围内的实数,决定其探索最优解空间的能力.在迭代优化过程中,粒子的位置和速度更新公式如下[23,24]:

其中,ωk为惯性权重;c1和c2是加速度常数,通常设置为2;r1,r2,ρ 均为0—1 之间的随机数;和分别为第i个质点在第k时刻的速度和位置;为单个粒子历史最优值和全体粒子历史最优值;K为最大迭代次数;ωini为初始惯性权重;ωend为迭代至最大次数时的惯性权重.同时,选取合适的粒子群规模(设置为5)和学习因子(设置为2)等算法参数,让该算法在优化的前期能更快地收敛并获得全局最优解.基于BPSO 算法的数字化仿真建模迭代寻优流程如图2(b)所示,两组算法的仿真模型及目标函数等参数设置均保持一致,唯一不同点在于两种算法迭代寻优的机制上.

3 仿真优化结果与讨论

上述两种算法的数字化仿真建模及计算均是在CPU 处理器为Intel Core i5-10400F@2.9 GHz,64G 内存,64 位Windows 10 操作系统下,通过软件Matlab 2018a 和FDTD 2020 互关联运行完成计算.DBS 和BPSO 两种算法优化后的单元结构和相应透射率曲线如图3 所示,图3(a)—(d)中两种算法的最终优化结构完全不同,是因为FOM 函数是一组非凸函数,具有多个局部最优点,两组算法均陷入了局部最优解.DBS 算法本身是局部优化算法,容易陷入局部最优解,而BPSO 算法是全局优化算法,在本工作中也陷入了局部最优解,可能是算法参数设置不合适导致粒子群收敛于局部最优值.对比图3(e)—(h)的透射率曲线,DBS 算法设计的插入损耗分别为–0.0582 dB,–0.0268 dB;而BPSO 算法设计的插入损耗分别为–0.0325 dB,–0.1070 dB.综合比较两种算法设计的器件性能,BPSO 算法寻优的解其带外抑制效果好、半峰全宽窄,但中心频率处透射率略低.其中插入损耗IL 的计算公式为 10 lgT,T为中心频率的透射率.

图3 DBS 和BPSO 两种算法优化后的单元结构和相应透射率曲线Fig.3.Unit structure and corresponding transmittance curve optimized by DBS and BPSO algorithms,respectively.

两种算法的FOM 值随迭代次数变化曲线如图4 所示,DBS 算法共需要迭代7 次,但单次迭代需要遍历全局共55 个像素点,故完成一次优化需要进行385 次FDTD 的仿真计算.而BPSO 算法共需要迭代65 次,但单次迭代需要计算粒子群中5 种粒子对应的不同结构,故完成一次优化需要进行325 次FDTD 的仿真计算.两种算法在FDTD仿真过程中均设置相同的尺寸大小、像素数和网格密度等参数,其单次仿真平均计算时间约为30 s,故DBS 和BPSO 算法分别完成一次优化的仿真时间为11550 s 和9750 s,即BPSO 算法计算时间短、优化效率高.

图4 (a) DBS 和(b) BPSO 两种算法在优化过程中FOM 值随迭代次数的变化Fig.4.Variation curve of FOM value versus iteration number in the optimization process of (a) DBS and (b) BPSO algorithms,respectively.

图4(a)中DBS 算法的迭代过程比较平缓且FOM 值最终收敛于0.069,而图4(b)中BPSO 算法的迭代过程首先快速下降,表现为BPSO 算法仅需60 次仿真即可让FOM 值下降至0.15,而DBS算法却需要220 次仿真,随后迭代次数增加下降趋势平缓,直到FOM 值几乎不发生变化,最终收敛于0.056.综上可见,BPSO 算法依据设定的目标功能和约束条件在整个解空间中进行迭代寻优,在全局优化过程中利用每组结构的历史最优值与全局最优值作为下一次的迭代方向,相较于DBS 算法在优化前期能更快地收敛.同时BPSO 算法通过不断调节优化区域的结构使其功能不断地逼近极值,最终收敛的FOM 值更小,因此该算法最终优化结果与预先设定的目标函数匹配程度高.

BPSO 算法通过优质的进化策略迭代寻优,用目标导向代替随机迭代,选取优化过程中不同FOM值观察算法的寻优机制.如图5(a)所示,当FOM值为0.150 时,结构的中心频率处于0.64 THz,半峰全宽为133.4 GHz,带外透射率幅值较高;经过数次稳定且有准确进化方向的迭代,FOM 值降低至0.100,其中心频率也移动至0.49 THz,半峰全宽变窄至35.6 GHz,但带外透射率幅值依旧较高;当FOM 值收敛于0.056 时,中心频率和半峰全宽逼近设置参数,且带外透射率幅值下降至0.05 以下.综上寻优过程,BPSO 算法首先满足中心频率及带宽的要求,随后抑制带外透射率以实现目标功能.为测试该智能化THz BPFs 逆向设计框架的稳定性,输入其他不同目标频率至该框架中,其优化结果如图5(b)所示,结果表明本文搭建的逆向设计框架可实现不同功能器件结构的精准求解.上述规范化设计流程不限于此案例,还适用于带阻、高通、低通、偏振等其他功能性太赫兹器件的设计,具有通用性.

图5 (a) BPSO 算法在迭代寻优过程中不同FOM 值对应的透射率曲线;(b) BPSO 算法设计不同中心频率处THz BPFs 的透射率曲线Fig.5.(a) Transmittance curves corresponding to different FOM values in optimization process of BPSO algorithm;(b) transmittance curves of THz BPFs at different center frequencies designed by BPSO algorithm.

图6(a),(b)分别为逆向设计和传统设计按目标需求设计得到的性能最优的带通滤波器,其对应的透射率光谱如图6(c),(d)所示.对比光谱结果可以看出逆向设计得到的结果fgoal处透射率更高、BFWHM更窄且带外抑制效果更好,更符合目标性能.

图6 逆向设计和传统设计优化后的单元结构及其透射率曲线Fig.6.Unit structure and its transmittance curve optimized by reverse design and traditional design.

逆向设计直接从目标出发,设计出的器件性能更加优异.相比于传统设计,通过逆向设计的精细结构变化对透射光谱有明显的影响.以图5(b)中,中心频率为0.6 THz 的带通滤波器为例,其结构的具体简化步骤如图7 所示,首先将图7(a)中悬空的金属去除(箭头所指),随后将图7(b)中锯齿状边缘用空气沿直线近似填补,最后去除图7(c)中红色圆圈所指的空气矩形,得到图7(d)所示的简化单元结构.从图7(e)—(h)可以看出,随着单元结构中几个像素单元的变化,其透射光谱发生着显著的变化.其中心频率从0.59 THz 逐渐移动 至0.85 THz,BFWHM从78.8 GHz 逐渐展宽至145.2 GHz,由此可以看出精细结构对带宽、透射率以及带外抑制具有深刻的影响,体现了逆向设计的优越性.值得注意的是,简化后中心频率处峰值透射率略微降低,带外透射率变高,同时简化后的样品透射光谱与传统十字型、Y 型、圆孔型FSS 的透射光谱相似.

图7 (a)—(d)结构简化过程;(e)—(h)相应的透射率光谱Fig.7.(a)–(d) Process of structural simplification;(e)–(h) corresponding transmittance spectra.

4 结论

本文构建了基于DBS 和BPSO 算法的逆向设计框架并以此设计了窄带、低插入损耗、带外抑制强的数字型THz BPFS.研究发现BPSO 算法相较于DBS 算法收敛速度更快、寻优能力更强、优化结果性能更优.通用化和标准化的逆向设计框架优于传统人工正向设计,可以满足当下多样化、快节奏的设计需求,有助于高性能和新奇功能器件的设计与优化.