声子晶体梁边界态的强化学习设计方案与本征模式分析

金亚斌，何良书

（同济大学航空航天与力学学院，上海 200092）

凝聚态物理中拓扑绝缘体的概念为声波/弹性波超构材料的研究注入了全新的活力，其典型特征是材料内部表现为绝缘状态，而表面具有受拓扑保护的边界态［1-2］。拓扑绝缘体可以为系统提供方向选择性传输［3］、非对称传输［4］、拓扑法诺共振［5］和高鲁棒［6］等新颖波动调控。拓扑保护边界态的实现通常依赖于系统带隙边缘的能带翻转引起上下模式交换进而产生拓扑相位。从普通相位到拓扑相位的过程无法通过连续变化完成转换，其间必然经历了带隙闭合形成简并点的状态［7-8］。对于一维系统，拓扑相位用Zak相位来表示，研究表明系统的反射相位与Zak相位之间有明确的关系，因而为了便于计算，常采用反射相位表征拓扑性质［9］。当一维声子晶体中间界面两侧的半结构具有不同拓扑性质时，在其界面处存在由能带翻转形成的边界态。然而，边界态频率的确定与半结构的反射系数、反射相位以及整体结构的本征模式之间的关系尚未明确，有待进一步探讨。

作为一维弹性波系统的代表之一，声子晶体梁结构是工程中重要的基本结构，其隔振［10］和减振［11］等功能被深入地研究，若结合拓扑态设计能够更好地实现弹性波的精准和奇异调控［12］。一般而言，边界态的鲁棒性与带隙的宽度相关，宽带隙结构形成的边界态在局域化程度和抗干扰能力等方面优于窄带隙结构［13］。因此，根据带隙范围设计结构以实现目标范围内的高质量弹性波传输对声子晶体梁的实际应用有重要意义。近年来，人工智能与各学科的交叉融合在不断加深，其子领域机器学习已在指导超构材料设计方面占据了重要地位［14］。在机器学习设计方案中可分为两种类型，一类较常见的方案是以数据驱动方式采用神经网络等机器学习方法［15-16］，实现目标响应与结构之间的非线性映射。另一类则是最近兴起的基于环境交互的强化学习方案［17-19］，通过刺激结构参数在参数空间中探索以获得具有预定目标响应的结构。后者不依赖于前期的数据获取，通过合理的分配智能体，设计探索策略来逐步寻求最优解。从应用角度而言，通常需要针对特定目标带隙进行结构设计。在此类设计研究中，Luo等［17］采用强化学习探索结构厚度设计了层状声子晶体；Wu等［18］则通过探索质量分布设计了一维弹性棒。上述研究充分验证了强化学习在超结构逆设计中的有效性，但尚未进一步探索所设计的结构除带隙特征外的其他性质，如拓扑相位。由于带隙的拓扑相位特征与边界态品质相关联且较为敏感，因此讨论强化学习所设计的结构是否能构造高品质的边界态具有一定探索意义。

本文研究工程中应用广泛的声子晶体梁结构，采用传递矩阵法推导系统的透反射谱和色散关系。在此基础上，借助强化学习方案以实现具有目标带隙范围的梁结构设计，并通过拓扑相变条件获得具有不同拓扑性质的两个半结构。进一步，通过本征模式计算条件获得整体结构的边界态频率，并讨论半结构的反射系数和反射相位与边界态频率之间的关系，从而验证本征模式条件计算的正确性。最后，讨论了所设计边界态的鲁棒性。

1 声子晶体梁模型分析

本文通过在铝制直梁上周期性地排列厚质量块，构造了如图1a所示的周期性阶梯截面声子晶体梁模型，其单胞包含两个相同的厚质量块，薄块和厚块交替排列，晶格常数a=2a1+2a2+a3，厚度h1=10 mm，h2=20 mm，宽度b=20 mm，其材料参数为：密度ρ=2 700 kg·m-3，杨氏模量E=70 GPa，泊松比ν=0.33。

考虑纵波在弹性介质中传播的控制方程，有：

图1 声子晶体梁结构以及有限系统透反射示意图Fig.1 Schematic diagram of phononic beam and finite system for transmission and reflection

式中：q=ω/cL为波数；A+和A-分别表示右行波和左行波波幅；ω为圆频率；i为虚数单位。在梁横截面面积突变处根据位移和力的连续条件，有：

式中：xj表示第j个界面的横坐标，当j是奇数（偶数）时，α=1（α=2），β=2（β=1）；S1和S2分别表示薄块和厚块的横截面积。对整体结构应用传递矩阵法，得：

另一方面，考虑计算色散关系时，在单胞中应用传递矩阵法，并在单胞两侧应用布洛赫定理形成周期性的无限系统条件。最终，求解整体传递矩阵与相因子满足的系数行列式可得色散关系。具体推导过程可参考文献［19］，其色散关系为

式中：k为Bloch波数。

2 强化学习逆设计具有预期带隙结构

2.1 强化学习模型建立

强化学习模型通常包括三个部分：智能体、环境和策略［20］。智能体是强化学习中的动作实体，处于某一环境中。智能体根据当前状态执行动作，从环境中获得与目标差异的反馈以评价该动作执行所获得的收益，并使用反馈更新策略，再进入下一状态。根据策略中储存的学习经验，可以迫使智能体在后续探索中遇到同种状态时选择合适的动作。Q学习［21］是强化学习的重要代表性算法之一，其状态价值函数称为Q（s，a）函数，储存在由状态s和动作a构成的二维Q表中。本文采用Q学习的设计方案逆设计具有预期带隙的结构，流程图如图2所示。

在图2中，将整个设计流程分为强化学习部分和环境部分。在强化学习部分，将单胞的三个长度参数视为智能体，其初始值为a1=a2=a3=30 mm，其余结构和几何参数不变。为保证探索过程处在一个合理范围内，此处给定参数限制范围为

a1∈[1，60]mm，a2∈[1，60]mm，a3∈[1，160]mm，设定每一步的增量d=1 mm。在每个状态下包含6个可选动作，即dA1：a'1=a1+d；dA2：a'1=a1-

d；dB1：a'2=a2+d；dB2：a'2=a2-d；dC1：a'3=a3+d；dC2：a'3=a3-d。因此，整个参数空间包含60×60×160=576 000种状态，若执行动作后，跳出参数范围，即不存在此种状态，则该动作无效。在开始阶段，设定目标带隙频率［m，n］（本文讨论第三带隙），并初始化Q表，表中所有值均为零。进入外层循环，对于每一个片段（episode），当某一状态下的智能体对应的目标带隙误差评估值Rerr1小于或等于给定的阈值时，结束当前片段，并判断是否达到了设定的最大训练轮次。若达到最大训练轮次，则训练结束，否则进入下一个片段。值得注意的是，若Rerr1=0，则表明所找到的结构对应的带隙与目标完全吻合。但就实际而言，上述长度参数的变化是按一定增量变化的，可能导致部分目标无法达成，且本征方程的求解存在一定舍入误差等因素，因此完全吻合条件过于理想化。需通过设定合理的阈值，使得智能体的探索在可接受范围内满足目标，训练方能顺利进行。对于某一片段的每一步（step），即内层循环，采用ε贪心策略选择动作实现探索和优化之间的平衡（ε介于0和1之间），以ε概率随机选择动作，代表了算法的探索性质，（1-ε）概率选择当前状态下Q（s，a）最大函数值对应的动作，代表了算法的优化性质。执行完动作之后，进入状态s'，以状态s'下所有动作的价值函数极大值来构造更新项（式（6））。为了增加训练的稳定性，本文采用梯度下降的贪心策略，初始ε设定为0.1，随训练进行，每5个片段ε变为原来的0.9倍。动作执行后，智能体将由a1，a2，a3变为a'1，a'2，a'3，二者对应的目标带隙误差评估值进行比较。若有所降低，则表明此前所选动作有助于接近目标，根据比较结果设定奖励值以更新Q表，此后进入下一状态，直至Rerr1≤阈值，内层循环跳出。图2中的两次更新Q表，根据Q学习的off-policy规则执行，即：

图2 设计具有预期带隙声子晶体梁的强化学习框架Fig.2 Reinforcement learning framework for design of phononic beam with anticipating bandgap

式中：α为学习率；γ为折减系数；r为奖励值；s和s'分别为执行动作前后的状态。

训练流程中的目标带隙误差评估值由环境部分提供，环境部分接收智能体参数，根据式（5）计算能带以获得当前的带隙频率［x，y］，则误差评估值可以用公式的形式给出。

2.2 强化学习模型验证

为验证2.1节中强化学习模型框架的效果，本节通过随机给定4个探索目标进行模型训练和测试。分别是第三带隙目标频率为［35 000，42 000］Hz，［40 000，48 000］Hz，［46 000，52 000］Hz和［50 000，60 000］Hz，设定阈值为100，最大训练片段数为35，学习率为0.01，折减系数为0.9。从图3中可以看出，对于上述4个目标的训练，每个训练片段达到目标所使用的步数都随着训练的进行逐渐趋于稳定的较小值。根据图2中的流程图，每一个片段最后都能探索获得满足阈值的结构，否则无法进入下一个片段，因此一个训练片段也是足够的。设置多个片段目的是探究模型逐渐找到一个稳定演化路径的过程，且所获得演化路径是较为节省探索步数的路径。

图3 每个训练片段（episode）达到目标所使用的步数（step）Fig.3 Number of steps used in each episode to achieve the goal

图4依次显示了4个探索目标最后一个训练片段的误差评估值变化图。从图4中可以看出，所有训练的误差在最后一个片段都呈波动下降，最终将低于阈值，其最终误差如表1所示。由于流程中为减小陷入局部极小的可能性，动作选择具有一定的随机性，且由于目标不唯一，即所有低于阈值的结果都是满足的目标。因此训练过程的误差不可避免存在一定波动性，但只需最终收敛至阈值以下，则表明满足目标的结构已探索成功。表1补充展现了训练的结果，所获得的结构参数对应的带隙均与原目标带隙十分接近。

表1 强化学习模型的训练结果Tab.1 Training results of reinforcement learning model

图4 最后一个训练片段（episode）中每一步（step）的误差评估值Fig.4 Error evaluation value of each step in the last episode

3 目标边界态设计与分析

在第2节中，对所提出的强化学习设计方案的有效性和稳健性进行了多个目标以及多个训练片段的重复测试。在此基础上，本节将进行边界态的设计与分析。由图2可知，阈值的大小控制着设计误差与探索时间的平衡。阈值越小，表示迭代退出条件越严格，通常所需探索时间将越长，但可能获得设计误差更小的结构。在本节中，由于仅需要获得目标结构，因此只进行一个训练片段即可退出。为使逆设计获得的结构更加精确，在本节中加强迭代退出条件，将阈值设置为10。给定预期目标带隙为［32 000，42 000］Hz，则阈值占所取带隙宽度的百分比仅为0.1%，因此在该阈值限制下能得到更为精确的设计。通过采用强化学习设计方案获得了满足目标的结构，相应的结构参数为a1=53 mm，a2=29 mm和a3=42 mm，则晶格常数a=206 mm。该组参数对应的第三带隙为［31 999，42 001］Hz，与预期目标几乎一致，误差评估值Rerr=1.4，表明了设计的精确性。

保证晶格常数和厚块的长度a2不变，改变a3的长度，即将两个厚块之间的距离由小到大变化，研究第三带隙边缘频率，如图5所示。随着a3增大，带隙经历了开-闭-开的过程，在a3=(a-2a2)/2=74 mm处（见中间竖直虚线）出现了拓扑过渡点，其左右两侧能带出现了上下翻转，虽然左右两侧能带可以完全对称，但模式的转换表明体系经历了拓扑相变，其左右两侧拓扑性质不同。利用其对称性，此前强化学习所获得结构A（左侧竖直虚线）存在与之带隙完全一致但拓扑性质不同的结构B（右侧竖直虚线），该结构对应的参数为a1=21 mm，a2=29 mm和a3=106 mm。将结构A与结构B拼接成整体的有限混合结构AB（如图6中的插图所示），结构A和B分别包含5个单胞，则当采用边界态频率激发时，能量将局域于结构的中间界面处。为了确定边界态频率，通常做法是将整体结构视为超元胞并计算其能带，带隙中的平带对应频率即为边界态频率。此方法中平带往往对应着一个频率范围，因此无法精确地获得边界态频率。

图5 第三带隙边缘频率与a3的关系Fig.5 3rd band gap edge frequency versus a3

图6 有限混合结构AB的传递矩阵中T11元素（abs（T11））与频率的关系Fig.6 T11 element of transfer matrix(abs(T11))of finite hybrid structure AB versus frequency

本文讨论从分析本征模式的角度一次性精准确定所有边界态频率。本征模式的一个显著特征是此时系统没有任何外部输入，该模式的本质是系统的固有状态，而边界态属于一种系统的固有状态。因此，通过对式（4）应用A+1=A-4N+1=0的边界条件，形成一个没有外部输入却存在输出的系统，得：

式（7）表明，仅通过关注系统传递矩阵的T11元素即可找到边界态。采用第1节的公式推导，可以建立系统的传递矩阵，并将T11的模与频率的关系绘制在图6中，其中灰色区域为带隙。可以看出，在第1和第3带隙中有明显的T11趋于零的频率出现，分别是11 660和36 680 Hz，均为边界态频率。第2和第4带隙中则没有边界态出现，因为其中不会出现拓扑相变，其左右两侧拓扑性质是相同的，这与文献［19］中的讨论分析是一致的。值得注意的是，从本征模式的角度分析可以同时找出多个边界态频率，若将频率继续提高，还会出现更高频的边界态频率。此外，第3带隙的宽度远大于第1带隙，其向下的尖峰更加尖锐，表明该边界态局域化程度和鲁棒性都将更高。

为验证本征模式条件所获边界态频率的准确性，根据第1节中的推导，计算有限混合结构AB的反射系数与透射系数，如图7a和7b所示，带隙中的峰和谷对应的频率与图6中的结果是一致的。图7c和7d分别展现了结构AB在11 660和36 680 Hz时的位移图，能量局域在整体结构的中部界面处。由于存在A+1=1的能量输入，结构的左侧边界存在一定范围的入射引起的波动。通过对比可以看出，图7c的体结构在远离中间界面时也明显存在较大的位移，表明第一带隙内的边界态局域化程度低于第三带隙，验证了根据图6的推测。

图7 结构AB的反射系数，透射系数和位移Fig.7 Reflection coefficient,transmission coefficient and displacement of structure AB

更进一步，分别研究半结构A和半结构B。图8a和图8b分别为半结构A和B的反射系数和反射相位。将两种结构的反射相位相加，可得到图8c。在边界态频率处，由于波的相消干涉，在界面处形成局域态，其反射相位之和为零，即prNA+prNB=0。因此，可假设此时反射系数分别为rNA=a+bi，rNB=a-bi，则rNA·rNB=a2+b2。由于带隙内反射系数的模接近于1，从而可得另一个确定边界态频率的条件rNA·rNB=1。在图8d中，36 680 Hz处的反射系数乘积比11 660 Hz处更接近于1，边界态条件得到了更严格的满足。因此，正如图7c和7d所展现的，36 680 Hz处的能量局域效果显著好于11 660 Hz处。以上两个条件表明不需要研究整体结构，仅研究两个半结构单独的反射情况也可以确定边界态频率。同样地，图8c和8d中三角形和圆形标记点的频率与本征模式所获边界态频率一致，从半结构的角度再次验证了本征模式所获边界态频率的准确性。

图8 半结构A和B的反射系数和反射相位Fig.8 Reflection coefficient and reflection phase of semi structures A and B

最后，以第3带隙形成的边界态为例研究其鲁棒性。考虑两种扰动：①对每个单胞的厚块沿x方向的位置进行随机位移，同一个单胞内两个厚块对称相互靠近或远离；②随机将若干（M个）厚块变成薄块形成缺陷。使用rNA·rNB=1的条件考量鲁棒性的强弱。对于扰动①，定义随机变量δ来量化位置的无序强度，其含义是各厚块的随机移动范围为［-δ，δ］。从图9a中可以看出，虽然位置的无序会引起一定的波动，但边界态仍然可以稳定地保持。在图9b中，随机缺陷同样无法有效地破坏边界态。以上两种检验结果表明，所构造的边界态具有较高的鲁棒性。对于纵波在对称梁中的传播，有两个对称平面xoy和xoz［4］。扰动①和扰动②均未破坏两个对称平面的对称性，因此纵波在该结构中的传播得到了有效的对称性保护。

图9 以rNA·rNB=1条件检验边界态的鲁棒性Fig.9 Robustness of the edge state judged by condition rNA·rNB=1

4 结语

本文设计和分析了声子晶体梁结构的边界态。首先，采用传递矩阵法推导获得了周期性阶梯截面声子晶体梁的透反射和色散关系。借助理论结果建立并验证了一种强化学习模型框架，以实现预期带隙设计梁结构的功能。进一步，采用该模型和拓扑理论设计了一个具有预期边界态的有限混合梁结构，并从本征模式的角度计算其边界态频率。在此基础上，分别从整体结构的透反射和两个半结构单独的透反射两种途径验证了本征模式计算边界态频率的准确性。最后，通过对厚块的位置扰动和厚块的缺陷两方面检验了所设计边界态的鲁棒性。以上研究工作力求将声子晶体梁从目标波动功能出发的设计过程简单化和程序化，并提出从本征模式的角度来分析边界态，全面地理解边界态与结构本身的关系。鉴于声子晶体梁在工程中具有优异的波动调控功能，该研究为声子晶体梁的设计与应用提供了一定的帮助。

作者贡献声明：

何良书：理论推导，数值计算，论文撰写。

金亚斌：研究方法，学术指导，论文修改。